diff --git a/python/fbprophet/diagnostics.py b/python/fbprophet/diagnostics.py
index 2f47bd3..17d129e 100644
--- a/python/fbprophet/diagnostics.py
+++ b/python/fbprophet/diagnostics.py
@@ -184,6 +184,7 @@ def prophet_copy(m, cutoff=None):
         weekly_seasonality=False,
         daily_seasonality=False,
         holidays=m.holidays,
+        seasonality_mode=m.seasonality_mode,
         seasonality_prior_scale=m.seasonality_prior_scale,
         changepoint_prior_scale=m.changepoint_prior_scale,
         holidays_prior_scale=m.holidays_prior_scale,
diff --git a/python/fbprophet/forecaster.py b/python/fbprophet/forecaster.py
index 87429ec..abf9a37 100644
--- a/python/fbprophet/forecaster.py
+++ b/python/fbprophet/forecaster.py
@@ -70,6 +70,7 @@ class Prophet(object):
         lower_window=-2 will include 2 days prior to the date as holidays. Also
         optionally can have a column prior_scale specifying the prior scale for
         that holiday.
+    seasonality_mode: 'additive' (default) or 'multiplicative'.
     seasonality_prior_scale: Parameter modulating the strength of the
         seasonality model. Larger values allow the model to fit larger seasonal
         fluctuations, smaller values dampen the seasonality. Can be specified
@@ -100,6 +101,7 @@ class Prophet(object):
             weekly_seasonality='auto',
             daily_seasonality='auto',
             holidays=None,
+            seasonality_mode='additive',
             seasonality_prior_scale=10.0,
             holidays_prior_scale=10.0,
             changepoint_prior_scale=0.05,
@@ -132,6 +134,7 @@ class Prophet(object):
             holidays['ds'] = pd.to_datetime(holidays['ds'])
         self.holidays = holidays
 
+        self.seasonality_mode = seasonality_mode
         self.seasonality_prior_scale = float(seasonality_prior_scale)
         self.changepoint_prior_scale = float(changepoint_prior_scale)
         self.holidays_prior_scale = float(holidays_prior_scale)
@@ -173,6 +176,10 @@ class Prophet(object):
                     raise ValueError('Holiday upper_window should be >= 0')
             for h in self.holidays['holiday'].unique():
                 self.validate_column_name(h, check_holidays=False)
+        if self.seasonality_mode not in ['additive', 'multiplicative']:
+            raise ValueError(
+                "seasonality_mode must be 'additive' or 'multiplicative'"
+            )
 
     def validate_column_name(self, name, check_holidays=True,
                              check_seasonalities=True, check_regressors=True):
@@ -189,7 +196,8 @@ class Prophet(object):
             raise ValueError('Name cannot contain "_delim_"')
         reserved_names = [
             'trend', 'additive_terms', 'daily', 'weekly', 'yearly',
-            'holidays', 'zeros', 'extra_regressors', 'yhat'
+            'holidays', 'zeros', 'extra_regressors_additive',
+            'extra_regressors_multiplicative', 'yhat',
         ]
         rn_l = [n + '_lower' for n in reserved_names]
         rn_u = [n + '_upper' for n in reserved_names]
@@ -410,6 +418,7 @@ class Prophet(object):
         -------
         holiday_features: pd.DataFrame with a column for each holiday.
         prior_scale_list: List of prior scales for each holiday column.
+        holiday_names: List of names of holidays
         """
         # Holds columns of our future matrix.
         expanded_holidays = defaultdict(lambda: np.zeros(dates.shape[0]))
@@ -461,9 +470,11 @@ class Prophet(object):
             prior_scales[h.split('_delim_')[0]]
             for h in holiday_features.columns
         ]
-        return holiday_features, prior_scale_list
+        return holiday_features, prior_scale_list, list(prior_scales.keys())
 
-    def add_regressor(self, name, prior_scale=None, standardize='auto'):
+    def add_regressor(
+        self, name, prior_scale=None, standardize='auto', mode=None
+    ):
         """Add an additional regressor to be used for fitting and predicting.
 
         The dataframe passed to `fit` and `predict` will have a column with the
@@ -472,6 +483,10 @@ class Prophet(object):
         coefficient is given a prior with the specified scale parameter.
         Decreasing the prior scale will add additional regularization. If no
         prior scale is provided, self.holidays_prior_scale will be used.
+        Mode can be specified as either 'additive' or 'multiplicative'. If not
+        specified, self.seasonality_mode will be used. 'additive' means the
+        effect of the regressor will be added to the trend, 'multiplicative'
+        means it will multiply the trend.
 
         Parameters
         ----------
@@ -481,6 +496,8 @@ class Prophet(object):
         standardize: optional, specify whether this regressor will be
             standardized prior to fitting. Can be 'auto' (standardize if not
             binary), True, or False.
+        mode: optional, 'additive' or 'multiplicative'. Defaults to
+            self.seasonality_mode.
 
         Returns
         -------
@@ -492,16 +509,23 @@ class Prophet(object):
         self.validate_column_name(name, check_regressors=False)
         if prior_scale is None:
             prior_scale = float(self.holidays_prior_scale)
+        if mode is None:
+            mode = self.seasonality_mode
         assert prior_scale > 0
+        if mode not in ['additive', 'multiplicative']:
+            raise ValueError("mode must be 'additive' or 'multiplicative'")
         self.extra_regressors[name] = {
             'prior_scale': prior_scale,
             'standardize': standardize,
             'mu': 0.,
             'std': 1.,
+            'mode': mode,
         }
         return self
 
-    def add_seasonality(self, name, period, fourier_order, prior_scale=None):
+    def add_seasonality(
+        self, name, period, fourier_order, prior_scale=None, mode=None
+    ):
         """Add a seasonal component with specified period, number of Fourier
         components, and prior scale.
 
@@ -514,12 +538,18 @@ class Prophet(object):
         seasonality_prior_scale provided on Prophet initialization (defaults
         to 10).
 
+        Mode can be specified as either 'additive' or 'multiplicative'. If not
+        specified, self.seasonality_mode will be used (defaults to additive).
+        Additive means the seasonality will be added to the trend,
+        multiplicative means it will multiply the trend.
+
         Parameters
         ----------
         name: string name of the seasonality component.
         period: float number of days in one period.
         fourier_order: int number of Fourier components to use.
-        prior_scale: float prior scale for this component.
+        prior_scale: optional float prior scale for this component.
+        mode: optional 'additive' or 'multiplicative'
 
         Returns
         -------
@@ -537,10 +567,15 @@ class Prophet(object):
             ps = float(prior_scale)
         if ps <= 0:
             raise ValueError('Prior scale must be > 0')
+        if mode is None:
+            mode = self.seasonality_mode
+        if mode not in ['additive', 'multiplicative']:
+            raise ValueError("mode must be 'additive' or 'multiplicative'")
         self.seasonalities[name] = {
             'period': period,
             'fourier_order': fourier_order,
             'prior_scale': ps,
+            'mode': mode,
         }
         return self
 
@@ -560,9 +595,12 @@ class Prophet(object):
         list of prior scales for each column of the features dataframe.
         Dataframe with indicators for which regression components correspond to
             which columns.
+        Dictionary with keys 'additive' and 'multiplicative' listing the
+            component names for each mode of seasonality.
         """
         seasonal_features = []
         prior_scales = []
+        modes = {'additive': [], 'multiplicative': []}
 
         # Seasonality features
         for name, props in self.seasonalities.items():
@@ -575,26 +613,120 @@ class Prophet(object):
             seasonal_features.append(features)
             prior_scales.extend(
                 [props['prior_scale']] * features.shape[1])
+            modes[props['mode']].append(name)
 
         # Holiday features
         if self.holidays is not None:
-            features, holiday_priors = self.make_holiday_features(df['ds'])
+            features, holiday_priors, holiday_names = (
+                self.make_holiday_features(df['ds'])
+            )
             seasonal_features.append(features)
             prior_scales.extend(holiday_priors)
+            modes[self.seasonality_mode].extend(holiday_names)
 
         # Additional regressors
         for name, props in self.extra_regressors.items():
             seasonal_features.append(pd.DataFrame(df[name]))
             prior_scales.append(props['prior_scale'])
+            modes[props['mode']].append(name)
 
+        # Dummy to prevent empty X
         if len(seasonal_features) == 0:
             seasonal_features.append(
                 pd.DataFrame({'zeros': np.zeros(df.shape[0])}))
             prior_scales.append(1.)
 
         seasonal_features = pd.concat(seasonal_features, axis=1)
-        component_cols = self.regressor_column_matrix(seasonal_features)
-        return seasonal_features, prior_scales, component_cols
+        component_cols, modes = self.regressor_column_matrix(
+            seasonal_features, modes
+        )
+        return seasonal_features, prior_scales, component_cols, modes
+
+    def regressor_column_matrix(self, seasonal_features, modes):
+        """Dataframe indicating which columns of the feature matrix correspond
+        to which seasonality/regressor components.
+
+        Includes combination components, like 'additive_terms'. These
+        combination components will be added to the 'modes' input.
+
+        Parameters
+        ----------
+        seasonal_features: Constructed seasonal features dataframe
+        modes: Dictionary with keys 'additive' and 'multiplicative' listing the
+            component names for each mode of seasonality.
+
+        Returns
+        -------
+        component_cols: A binary indicator dataframe with columns seasonal
+            components and rows columns in seasonal_features. Entry is 1 if
+            that columns is used in that component.
+        modes: Updated input with combination components.
+        """
+        components = pd.DataFrame({
+            'col': np.arange(seasonal_features.shape[1]),
+            'component': [
+                x.split('_delim_')[0] for x in seasonal_features.columns
+            ],
+        })
+        # Add total for holidays
+        if self.holidays is not None:
+            components = self.add_group_component(
+                components, 'holidays', self.holidays['holiday'].unique())
+        # Add totals additive and multiplicative components, and regressors
+        for mode in ['additive', 'multiplicative']:
+            components = self.add_group_component(
+                components, mode + '_terms', modes[mode]
+            )
+            regressors_by_mode = [
+                r for r, props in self.extra_regressors.items()
+                if props['mode'] == mode
+            ]
+            components = self.add_group_component(
+                components, 'extra_regressors_' + mode, regressors_by_mode)
+            # Add combination components to modes
+            modes[mode].append(mode + '_terms')
+            modes[mode].append('extra_regressors_' + mode)
+        # Convert to a binary matrix
+        component_cols = pd.crosstab(
+            components['col'], components['component'],
+        )
+        # Add columns for additive and multiplicative terms, if missing
+        for name in ['additive_terms', 'multiplicative_terms']:
+            if name not in component_cols:
+                component_cols[name] = 0
+        # Remove the placeholder
+        component_cols.drop('zeros', axis=1, inplace=True, errors='ignore')
+        # Validation
+        if (
+            max(component_cols['additive_terms']
+            + component_cols['multiplicative_terms']) > 1
+        ):
+            raise Exception('A bug occurred in seasonal components.')
+        # Compare to the training, if set.
+        if self.train_component_cols is not None:
+            component_cols = component_cols[self.train_component_cols.columns]
+            if not component_cols.equals(self.train_component_cols):
+                raise Exception('A bug occurred in constructing regressors.')
+        return component_cols, modes
+
+    def add_group_component(self, components, name, group):
+        """Adds a component with given name that contains all of the components
+        in group.
+
+        Parameters
+        ----------
+        components: Dataframe with components.
+        name: Name of new group component.
+        group: List of components that form the group.
+
+        Returns
+        -------
+        Dataframe with components.
+        """
+        new_comp = components[components['component'].isin(set(group))].copy()
+        new_comp['component'] = name
+        components = components.append(new_comp)
+        return components
 
     def parse_seasonality_args(self, name, arg, auto_disable, default_order):
         """Get number of fourier components for built-in seasonalities.
@@ -656,6 +788,7 @@ class Prophet(object):
                 'period': 365.25,
                 'fourier_order': fourier_order,
                 'prior_scale': self.seasonality_prior_scale,
+                'mode': self.seasonality_mode,
             }
 
         # Weekly seasonality
@@ -668,6 +801,7 @@ class Prophet(object):
                 'period': 7,
                 'fourier_order': fourier_order,
                 'prior_scale': self.seasonality_prior_scale,
+                'mode': self.seasonality_mode,
             }
 
         # Daily seasonality
@@ -680,6 +814,7 @@ class Prophet(object):
                 'period': 1,
                 'fourier_order': fourier_order,
                 'prior_scale': self.seasonality_prior_scale,
+                'mode': self.seasonality_mode,
             }
 
     @staticmethod
@@ -785,7 +920,7 @@ class Prophet(object):
         history = self.setup_dataframe(history, initialize_scales=True)
         self.history = history
         self.set_auto_seasonalities()
-        seasonal_features, prior_scales, component_cols = (
+        seasonal_features, prior_scales, component_cols, _ = (
             self.make_all_seasonality_features(history))
         self.train_component_cols = component_cols
 
@@ -802,6 +937,8 @@ class Prophet(object):
             'sigmas': prior_scales,
             'tau': self.changepoint_prior_scale,
             'trend_indicator': int(self.growth == 'logistic'),
+            's_a': component_cols['additive_terms'],
+            's_m': component_cols['multiplicative_terms'],
         }
 
         if self.growth == 'linear':
@@ -891,7 +1028,10 @@ class Prophet(object):
             cols.append('floor')
         # Add in forecast components
         df2 = pd.concat((df[cols], intervals, seasonal_components), axis=1)
-        df2['yhat'] = df2['trend'] + df2['additive_terms']
+        df2['yhat'] = (
+            df2['trend'] * (1 + df2['multiplicative_terms'])
+            + df2['additive_terms']
+        )
         return df2
 
     @staticmethod
@@ -991,7 +1131,7 @@ class Prophet(object):
         -------
         Dataframe with seasonal components.
         """
-        seasonal_features, _, component_cols = (
+        seasonal_features, _, component_cols, modes = (
             self.make_all_seasonality_features(df)
         )
         lower_p = 100 * (1.0 - self.interval_width) / 2
@@ -1002,7 +1142,9 @@ class Prophet(object):
         for component in component_cols.columns:
             beta_c = self.params['beta'] * component_cols[component].values
 
-            comp = np.matmul(X, beta_c.transpose()) * self.y_scale
+            comp = np.matmul(X, beta_c.transpose())
+            if component in modes['additive']:
+                 comp *= self.y_scale
             data[component] = np.nanmean(comp, axis=1)
             data[component + '_lower'] = np.nanpercentile(
                 comp, lower_p, axis=1,
@@ -1012,54 +1154,6 @@ class Prophet(object):
             )
         return pd.DataFrame(data)
 
-    def regressor_column_matrix(self, seasonal_features):
-        components = pd.DataFrame({
-            'col': np.arange(seasonal_features.shape[1]),
-            'component': [x.split('_delim_')[0] for x in seasonal_features.columns],
-        })
-        # Add total for all additive components
-        components = components.append(pd.DataFrame({
-            'col': np.arange(seasonal_features.shape[1]),
-            'component': 'additive_terms',
-        }))
-        # Add totals for holidays and extra regressors
-        if self.holidays is not None:
-            components = self.add_group_component(
-                components, 'holidays', self.holidays['holiday'].unique())
-        components = self.add_group_component(
-            components, 'extra_regressors', self.extra_regressors.keys())
-        # Remove the placeholder
-        components = components[components['component'] != 'zeros']
-        # Convert to a binary matrix
-        component_cols = pd.crosstab(
-            components['col'], components['component'],
-        )
-        # Compare to the training, if set.
-        if self.train_component_cols is not None:
-            component_cols = component_cols[self.train_component_cols.columns]
-            if not component_cols.equals(self.train_component_cols):
-                raise Exception('A bug occurred in constructing regressors.')
-        return component_cols
-
-    def add_group_component(self, components, name, group):
-        """Adds a component with given name that contains all of the components
-        in group.
-
-        Parameters
-        ----------
-        components: Dataframe with components.
-        name: Name of new group component.
-        group: List of components that form the group.
-
-        Returns
-        -------
-        Dataframe with components.
-        """
-        new_comp = components[components['component'].isin(set(group))].copy()
-        new_comp['component'] = name
-        components = components.append(new_comp)
-        return components
-
     def sample_posterior_predictive(self, df):
         """Prophet posterior predictive samples.
 
@@ -1069,7 +1163,8 @@ class Prophet(object):
 
         Returns
         -------
-        Dictionary with posterior predictive samples for each component.
+        Dictionary with posterior predictive samples for the forecast yhat and
+        for the trend component.
         """
         n_iterations = self.params['k'].shape[0]
         samp_per_iter = max(1, int(np.ceil(
@@ -1077,12 +1172,20 @@ class Prophet(object):
         )))
 
         # Generate seasonality features once so we can re-use them.
-        seasonal_features, _, _ = self.make_all_seasonality_features(df)
+        seasonal_features, _, component_cols, _ = (
+            self.make_all_seasonality_features(df)
+        )
 
-        sim_values = {'yhat': [], 'trend': [], 'seasonal': []}
+        sim_values = {'yhat': [], 'trend': []}
         for i in range(n_iterations):
             for _j in range(samp_per_iter):
-                sim = self.sample_model(df, seasonal_features, i)
+                sim = self.sample_model(
+                    df=df,
+                    seasonal_features=seasonal_features,
+                    iteration=i,
+                    s_a=component_cols['additive_terms'],
+                    s_m=component_cols['multiplicative_terms'],
+                )
                 for key in sim_values:
                     sim_values[key].append(sim[key])
         for k, v in sim_values.items():
@@ -1099,9 +1202,8 @@ class Prophet(object):
 
         Returns
         -------
-        Dictionary with keys "trend", "seasonal", and "yhat" containing
-        posterior predictive samples for that component. "seasonal" is the sum
-        of seasonalities, holidays, and added regressors.
+        Dictionary with keys "trend" and "yhat" containing
+        posterior predictive samples for that component.
         """
         df = self.setup_dataframe(df.copy())
         sim_values = self.sample_posterior_predictive(df)
@@ -1132,7 +1234,7 @@ class Prophet(object):
 
         return pd.DataFrame(series)
 
-    def sample_model(self, df, seasonal_features, iteration):
+    def sample_model(self, df, seasonal_features, iteration, s_a, s_m):
         """Simulate observations from the extrapolated generative model.
 
         Parameters
@@ -1143,20 +1245,22 @@ class Prophet(object):
 
         Returns
         -------
-        Dataframe with trend, seasonality, and yhat, each like df['t'].
+        Dataframe with trend and yhat, each like df['t'].
         """
         trend = self.sample_predictive_trend(df, iteration)
 
         beta = self.params['beta'][iteration]
-        seasonal = np.matmul(seasonal_features.as_matrix(), beta) * self.y_scale
+        Xb_a = (
+            np.matmul(seasonal_features.as_matrix(), beta * s_a) * self.y_scale
+        )
+        Xb_m = np.matmul(seasonal_features.as_matrix(), beta * s_m)
 
         sigma = self.params['sigma_obs'][iteration]
         noise = np.random.normal(0, sigma, df.shape[0]) * self.y_scale
 
         return pd.DataFrame({
-            'yhat': trend + seasonal + noise,
-            'trend': trend,
-            'seasonal': seasonal,
+            'yhat': trend * (1 + Xb_m) + Xb_a + noise,
+            'trend': trend
         })
 
     def sample_predictive_trend(self, df, iteration):
diff --git a/python/fbprophet/plot.py b/python/fbprophet/plot.py
index a87222a..e98fc3d 100644
--- a/python/fbprophet/plot.py
+++ b/python/fbprophet/plot.py
@@ -77,7 +77,8 @@ def plot_components(
     """Plot the Prophet forecast components.
 
     Will plot whichever are available of: trend, holidays, weekly
-    seasonality, and yearly seasonality.
+    seasonality, yearly seasonality, and additive and multiplicative extra
+    regressors.
 
     Parameters
     ----------
@@ -103,8 +104,12 @@ def plot_components(
         components.append('holidays')
     components.extend([name for name in m.seasonalities
                     if name in fcst])
-    if len(m.extra_regressors) > 0 and 'extra_regressors' in fcst:
-        components.append('extra_regressors')
+    regressors = {'additive': False, 'multiplicative': False}
+    for name, props in m.extra_regressors.items():
+        regressors[props['mode']] = True
+    for mode in ['additive', 'multiplicative']:
+        if regressors[mode] and 'extra_regressors_{}'.format(mode) in fcst:
+            components.append('extra_regressors_{}'.format(mode))
     npanel = len(components)
 
     fig, axes = plt.subplots(npanel, 1, facecolor='w',
@@ -119,11 +124,6 @@ def plot_components(
                 m=m, fcst=fcst, name='trend', ax=ax, uncertainty=uncertainty,
                 plot_cap=plot_cap,
             )
-        elif plot_name == 'holidays':
-            plot_forecast_component(
-                m=m, fcst=fcst, name='holidays', ax=ax,
-                uncertainty=uncertainty, plot_cap=False,
-            )
         elif plot_name == 'weekly':
             plot_weekly(
                 m=m, ax=ax, uncertainty=uncertainty, weekly_start=weekly_start,
@@ -132,10 +132,14 @@ def plot_components(
             plot_yearly(
                 m=m, ax=ax, uncertainty=uncertainty, yearly_start=yearly_start,
             )
-        elif plot_name == 'extra_regressors':
+        elif plot_name in [
+            'holidays',
+            'extra_regressors_additive',
+            'extra_regressors_multiplicative',
+        ]:
             plot_forecast_component(
-                m=m, fcst=fcst, name='extra_regressors', ax=ax,
-                uncertainty=uncertainty, plot_cap=False,
+                m=m, fcst=fcst, name=plot_name, ax=ax, uncertainty=uncertainty,
+                plot_cap=False,
             )
         else:
             plot_seasonality(
@@ -242,7 +246,7 @@ def plot_weekly(m, ax=None, uncertainty=True, weekly_start=0):
     ax.set_xticks(range(len(days)))
     ax.set_xticklabels(days)
     ax.set_xlabel('Day of week')
-    ax.set_ylabel('weekly')
+    ax.set_ylabel('weekly ({})'.format(m.seasonalities['weekly']['mode']))
     return artists
 
 
@@ -284,7 +288,7 @@ def plot_yearly(m, ax=None, uncertainty=True, yearly_start=0):
         lambda x, pos=None: '{dt:%B} {dt.day}'.format(dt=num2date(x))))
     ax.xaxis.set_major_locator(months)
     ax.set_xlabel('Day of year')
-    ax.set_ylabel('yearly')
+    ax.set_ylabel('yearly ({})'.format(m.seasonalities['yearly']['mode']))
     return artists
 
 
@@ -334,7 +338,7 @@ def plot_seasonality(m, name, ax=None, uncertainty=True):
     ax.xaxis.set_major_formatter(FuncFormatter(
         lambda x, pos=None: fmt_str.format(dt=num2date(x))))
     ax.set_xlabel('ds')
-    ax.set_ylabel(name)
+    ax.set_ylabel('{} ({})'.format(name, m.seasonalities[name]['mode']))
     return artists
 
 
diff --git a/python/fbprophet/tests/test_diagnostics.py b/python/fbprophet/tests/test_diagnostics.py
index 4a8c154..2ea3f80 100644
--- a/python/fbprophet/tests/test_diagnostics.py
+++ b/python/fbprophet/tests/test_diagnostics.py
@@ -184,6 +184,7 @@ class TestDiagnostics(TestCase):
             [True, False],  # weekly_seasonality
             [True, False],  # daily_seasonality
             [None, holiday],  # holidays
+            ['additive', 'multiplicative'],  # seasonality_mode
             [1.1],  # seasonality_prior_scale
             [1.1],  # holidays_prior_scale
             [0.1],  # changepoint_prior_scale
@@ -214,6 +215,7 @@ class TestDiagnostics(TestCase):
                 self.assertEqual(m1.holidays, m2.holidays)
             else:
                 self.assertTrue((m1.holidays == m2.holidays).values.all())
+            self.assertEqual(m1.seasonality_mode, m2.seasonality_mode)
             self.assertEqual(m1.seasonality_prior_scale, m2.seasonality_prior_scale)
             self.assertEqual(m1.changepoint_prior_scale, m2.changepoint_prior_scale)
             self.assertEqual(m1.holidays_prior_scale, m2.holidays_prior_scale)
diff --git a/python/fbprophet/tests/test_prophet.py b/python/fbprophet/tests/test_prophet.py
index 35a115e..993d6d5 100644
--- a/python/fbprophet/tests/test_prophet.py
+++ b/python/fbprophet/tests/test_prophet.py
@@ -261,11 +261,12 @@ class TestProphet(TestCase):
         df = pd.DataFrame({
             'ds': pd.date_range('2016-12-20', '2016-12-31')
         })
-        feats, priors = model.make_holiday_features(df['ds'])
+        feats, priors, names = model.make_holiday_features(df['ds'])
         # 11 columns generated even though only 8 overlap
         self.assertEqual(feats.shape, (df.shape[0], 2))
         self.assertEqual((feats.sum(0) - np.array([1.0, 1.0])).sum(), 0)
         self.assertEqual(priors, [10., 10.])  # Default prior
+        self.assertEqual(names, ['xmas'])
 
         holidays = pd.DataFrame({
             'ds': pd.to_datetime(['2016-12-25']),
@@ -273,10 +274,12 @@ class TestProphet(TestCase):
             'lower_window': [-1],
             'upper_window': [10],
         })
-        feats, priors = Prophet(holidays=holidays).make_holiday_features(df['ds'])
+        m = Prophet(holidays=holidays)
+        feats, priors, names = m.make_holiday_features(df['ds'])
         # 12 columns generated even though only 8 overlap
         self.assertEqual(feats.shape, (df.shape[0], 12))
         self.assertEqual(priors, list(10. * np.ones(12)))
+        self.assertEqual(names, ['xmas'])
         # Check prior specifications
         holidays = pd.DataFrame({
             'ds': pd.to_datetime(['2016-12-25', '2017-12-25']),
@@ -285,8 +288,10 @@ class TestProphet(TestCase):
             'upper_window': [0, 0],
             'prior_scale': [5., 5.],
         })
-        feats, priors = Prophet(holidays=holidays).make_holiday_features(df['ds'])
+        m = Prophet(holidays=holidays)
+        feats, priors, names = m.make_holiday_features(df['ds'])
         self.assertEqual(priors, [5., 5.])
+        self.assertEqual(names, ['xmas'])
         # 2 different priors
         holidays2 = pd.DataFrame({
             'ds': pd.to_datetime(['2012-06-06', '2013-06-06']),
@@ -296,8 +301,10 @@ class TestProphet(TestCase):
             'prior_scale': [8] * 2,
         })
         holidays2 = pd.concat((holidays, holidays2))
-        feats, priors = Prophet(holidays=holidays2).make_holiday_features(df['ds'])
+        m = Prophet(holidays=holidays2)
+        feats, priors, names = m.make_holiday_features(df['ds'])
         self.assertEqual(priors, [8., 8., 5., 5.])
+        self.assertEqual(set(names), {'xmas', 'seans-bday'})
         holidays2 = pd.DataFrame({
             'ds': pd.to_datetime(['2012-06-06', '2013-06-06']),
             'holiday': ['seans-bday'] * 2,
@@ -305,7 +312,7 @@ class TestProphet(TestCase):
             'upper_window': [1] * 2,
         })
         holidays2 = pd.concat((holidays, holidays2))
-        feats, priors = Prophet(
+        feats, priors, names = Prophet(
             holidays=holidays2, holidays_prior_scale=4
         ).make_holiday_features(df['ds'])
         self.assertEqual(priors, [4., 4., 5., 5.])
@@ -357,8 +364,15 @@ class TestProphet(TestCase):
         self.assertEqual(m.weekly_seasonality, 'auto')
         m.fit(train)
         self.assertIn('weekly', m.seasonalities)
-        self.assertEqual(m.seasonalities['weekly'],
-                         {'period': 7, 'fourier_order': 3, 'prior_scale': 10.})
+        self.assertEqual(
+            m.seasonalities['weekly'],
+            {
+                'period': 7,
+                'fourier_order': 3,
+                'prior_scale': 10.,
+                'mode': 'additive',
+            },
+        )
         # Should be disabled due to too short history
         N = 9
         train = DATA.head(N)
@@ -375,8 +389,15 @@ class TestProphet(TestCase):
         self.assertNotIn('weekly', m.seasonalities)
         m = Prophet(weekly_seasonality=2, seasonality_prior_scale=3.)
         m.fit(DATA)
-        self.assertEqual(m.seasonalities['weekly'],
-                         {'period': 7, 'fourier_order': 2, 'prior_scale': 3.})
+        self.assertEqual(
+            m.seasonalities['weekly'],
+            {
+                'period': 7,
+                'fourier_order': 2,
+                'prior_scale': 3.,
+                'mode': 'additive',
+            },
+        )
 
     def test_auto_yearly_seasonality(self):
         # Should be enabled
@@ -386,7 +407,12 @@ class TestProphet(TestCase):
         self.assertIn('yearly', m.seasonalities)
         self.assertEqual(
             m.seasonalities['yearly'],
-            {'period': 365.25, 'fourier_order': 10, 'prior_scale': 10.},
+            {
+                'period': 365.25,
+                'fourier_order': 10,
+                'prior_scale': 10.,
+                'mode': 'additive',
+            },
         )
         # Should be disabled due to too short history
         N = 240
@@ -401,7 +427,12 @@ class TestProphet(TestCase):
         m.fit(DATA)
         self.assertEqual(
             m.seasonalities['yearly'],
-            {'period': 365.25, 'fourier_order': 7, 'prior_scale': 3.},
+            {
+                'period': 365.25,
+                'fourier_order': 7,
+                'prior_scale': 3.,
+                'mode': 'additive',
+            },
         )
 
     def test_auto_daily_seasonality(self):
@@ -410,8 +441,15 @@ class TestProphet(TestCase):
         self.assertEqual(m.daily_seasonality, 'auto')
         m.fit(DATA2)
         self.assertIn('daily', m.seasonalities)
-        self.assertEqual(m.seasonalities['daily'],
-                         {'period': 1, 'fourier_order': 4, 'prior_scale': 10.})
+        self.assertEqual(
+            m.seasonalities['daily'],
+            {
+                'period': 1,
+                'fourier_order': 4,
+                'prior_scale': 10.,
+                'mode': 'additive',
+            },
+        )
         # Should be disabled due to too short history
         N = 430
         train = DATA2.head(N)
@@ -423,8 +461,15 @@ class TestProphet(TestCase):
         self.assertIn('daily', m.seasonalities)
         m = Prophet(daily_seasonality=7, seasonality_prior_scale=3.)
         m.fit(DATA2)
-        self.assertEqual(m.seasonalities['daily'],
-                         {'period': 1, 'fourier_order': 7, 'prior_scale': 3.})
+        self.assertEqual(
+            m.seasonalities['daily'],
+            {
+                'period': 1,
+                'fourier_order': 7,
+                'prior_scale': 3.,
+                'mode': 'additive',
+            },
+        )
         m = Prophet()
         m.fit(DATA)
         self.assertNotIn('daily', m.seasonalities)
@@ -448,24 +493,38 @@ class TestProphet(TestCase):
         m = Prophet(holidays=holidays)
         m.add_seasonality(name='monthly', period=30, fourier_order=5,
                           prior_scale=2.)
-        self.assertEqual(m.seasonalities['monthly'],
-                         {'period': 30, 'fourier_order': 5, 'prior_scale': 2.})
+        self.assertEqual(
+            m.seasonalities['monthly'],
+            {
+                'period': 30,
+                'fourier_order': 5,
+                'prior_scale': 2.,
+                'mode': 'additive',
+            },
+        )
         with self.assertRaises(ValueError):
             m.add_seasonality(name='special_day', period=30, fourier_order=5)
         with self.assertRaises(ValueError):
             m.add_seasonality(name='trend', period=30, fourier_order=5)
         m.add_seasonality(name='weekly', period=30, fourier_order=5)
         # Test priors
-        m = Prophet(holidays=holidays, yearly_seasonality=False)
+        m = Prophet(
+            holidays=holidays, yearly_seasonality=False,
+            seasonality_mode='multiplicative',
+        )
         m.add_seasonality(name='monthly', period=30, fourier_order=5,
-                          prior_scale=2.)
+                          prior_scale=2., mode='additive')
         m.fit(DATA.copy())
-        seasonal_features, prior_scales, component_cols = (
+        self.assertEqual(m.seasonalities['monthly']['mode'], 'additive')
+        self.assertEqual(m.seasonalities['weekly']['mode'], 'multiplicative')
+        seasonal_features, prior_scales, component_cols, modes = (
             m.make_all_seasonality_features(m.history)
         )
         self.assertEqual(sum(component_cols['monthly']), 10)
         self.assertEqual(sum(component_cols['special_day']), 1)
         self.assertEqual(sum(component_cols['weekly']), 6)
+        self.assertEqual(sum(component_cols['additive_terms']), 10)
+        self.assertEqual(sum(component_cols['multiplicative_terms']), 7)
         if seasonal_features.columns[0] == 'monthly_delim_1':
             true = [2.] * 10 + [10.] * 6 + [4.]
             self.assertEqual(sum(component_cols['monthly'][:10]), 10)
@@ -480,10 +539,14 @@ class TestProphet(TestCase):
         m = Prophet()
         m.add_regressor('binary_feature', prior_scale=0.2)
         m.add_regressor('numeric_feature', prior_scale=0.5)
+        m.add_regressor(
+            'numeric_feature2', prior_scale=0.5, mode='multiplicative'
+        )
         m.add_regressor('binary_feature2', standardize=True)
         df = DATA.copy()
         df['binary_feature'] = [0] * 255 + [1] * 255
         df['numeric_feature'] = range(510)
+        df['numeric_feature2'] = range(510)
         with self.assertRaises(ValueError):
             # Require all regressors in df
             m.fit(df)
@@ -492,7 +555,13 @@ class TestProphet(TestCase):
         # Check that standardizations are correctly set
         self.assertEqual(
             m.extra_regressors['binary_feature'],
-            {'prior_scale': 0.2, 'mu': 0, 'std': 1, 'standardize': 'auto'},
+            {
+                'prior_scale': 0.2,
+                'mu': 0,
+                'std': 1,
+                'standardize': 'auto',
+                'mode': 'additive',
+            },
         )
         self.assertEqual(
             m.extra_regressors['numeric_feature']['prior_scale'], 0.5)
@@ -500,6 +569,8 @@ class TestProphet(TestCase):
             m.extra_regressors['numeric_feature']['mu'], 254.5)
         self.assertAlmostEqual(
             m.extra_regressors['numeric_feature']['std'], 147.368585, places=5)
+        self.assertEqual(
+            m.extra_regressors['numeric_feature2']['mode'], 'multiplicative')
         self.assertEqual(
             m.extra_regressors['binary_feature2']['prior_scale'], 10.)
         self.assertAlmostEqual(
@@ -512,10 +583,10 @@ class TestProphet(TestCase):
         self.assertAlmostEqual(df2['numeric_feature'][0], -1.726962, places=4)
         self.assertAlmostEqual(df2['binary_feature2'][0], 2.022859, places=4)
         # Check that feature matrix and prior scales are correctly constructed
-        seasonal_features, prior_scales, component_cols = (
+        seasonal_features, prior_scales, component_cols, modes = (
             m.make_all_seasonality_features(df2)
         )
-        self.assertEqual(seasonal_features.shape[1], 29)
+        self.assertEqual(seasonal_features.shape[1], 30)
         names = ['binary_feature', 'numeric_feature', 'binary_feature2']
         true_priors = [0.2, 0.5, 10.]
         for i, name in enumerate(names):
@@ -530,24 +601,35 @@ class TestProphet(TestCase):
             'ds': ['2014-06-01'],
             'binary_feature': [0],
             'numeric_feature': [10],
+            'numeric_feature2': [10],
         })
         with self.assertRaises(ValueError):
             m.predict(future)
         future['binary_feature2'] = 0
         fcst = m.predict(future)
-        self.assertEqual(fcst.shape[1], 28)
+        self.assertEqual(fcst.shape[1], 37)
         self.assertEqual(fcst['binary_feature'][0], 0)
         self.assertAlmostEqual(
-            fcst['extra_regressors'][0],
+            fcst['extra_regressors_additive'][0],
             fcst['numeric_feature'][0] + fcst['binary_feature2'][0],
         )
+        self.assertAlmostEqual(
+            fcst['extra_regressors_multiplicative'][0],
+            fcst['numeric_feature2'][0],
+        )
         self.assertAlmostEqual(
             fcst['additive_terms'][0],
-            fcst['yearly'][0] + fcst['weekly'][0] + fcst['extra_regressors'][0]
+            fcst['yearly'][0] + fcst['weekly'][0]
+                + fcst['extra_regressors_additive'][0],
+        )
+        self.assertAlmostEqual(
+            fcst['multiplicative_terms'][0],
+            fcst['extra_regressors_multiplicative'][0],
         )
         self.assertAlmostEqual(
             fcst['yhat'][0],
-            fcst['trend'][0] + fcst['additive_terms'][0],
+            fcst['trend'][0] * (1 + fcst['multiplicative_terms'][0])
+                + fcst['additive_terms'][0],
         )
         # Check fails if constant extra regressor
         df['constant_feature'] = 5
@@ -555,3 +637,46 @@ class TestProphet(TestCase):
         m.add_regressor('constant_feature')
         with self.assertRaises(ValueError):
             m.fit(df.copy())
+
+    def test_set_seasonality_mode(self):
+        # Setting attribute
+        m = Prophet()
+        self.assertEqual(m.seasonality_mode, 'additive')
+        m = Prophet(seasonality_mode='multiplicative')
+        self.assertEqual(m.seasonality_mode, 'multiplicative')
+        with self.assertRaises(ValueError):
+            Prophet(seasonality_mode='batman')
+
+    def test_seasonality_modes(self):
+        # Model with holidays, seasonalities, and extra regressors
+        holidays = pd.DataFrame({
+            'ds': pd.to_datetime(['2016-12-25']),
+            'holiday': ['xmas'],
+            'lower_window': [-1],
+            'upper_window': [0],
+        })
+        m = Prophet(seasonality_mode='multiplicative', holidays=holidays)
+        m.add_seasonality('monthly', period=30, mode='additive', fourier_order=3)
+        m.add_regressor('binary_feature', mode='additive')
+        m.add_regressor('numeric_feature')
+        # Construct seasonal features
+        df = DATA.copy()
+        df['binary_feature'] = [0] * 255 + [1] * 255
+        df['numeric_feature'] = range(510)
+        df = m.setup_dataframe(df, initialize_scales=True)
+        m.history = df.copy()
+        m.set_auto_seasonalities()
+        seasonal_features, prior_scales, component_cols, modes = (
+            m.make_all_seasonality_features(df))
+        self.assertEqual(sum(component_cols['additive_terms']), 7)
+        self.assertEqual(sum(component_cols['multiplicative_terms']), 29)
+        self.assertEqual(
+            set(modes['additive']),
+            {'monthly', 'binary_feature', 'additive_terms',
+             'extra_regressors_additive'},
+        )
+        self.assertEqual(
+            set(modes['multiplicative']),
+            {'weekly', 'yearly', 'xmas', 'numeric_feature',
+             'multiplicative_terms', 'extra_regressors_multiplicative'},
+        )
diff --git a/python/stan/unix/prophet.stan b/python/stan/unix/prophet.stan
index 746b0cb..4319b6f 100644
--- a/python/stan/unix/prophet.stan
+++ b/python/stan/unix/prophet.stan
@@ -85,6 +85,8 @@ data {
   vector[K] sigmas;     // Scale on seasonality prior
   real<lower=0> tau;    // Scale on changepoints prior
   int trend_indicator;  // 0 for linear, 1 for logistic
+  vector[K] s_a;        // Indicator of additive features
+  vector[K] s_m;        // Indicator of multiplicative features
 }
 
 transformed data {
@@ -102,12 +104,17 @@ parameters {
 
 transformed parameters {
   vector[T] trend;
+  vector[T] Xb_a;
+  vector[T] Xb_m;
 
   if (trend_indicator == 0) {
     trend = linear_trend(k, m, delta, t, A, t_change);
   } else if (trend_indicator == 1) {
     trend = logistic_trend(k, m, delta, t, cap, A, t_change, S);
   }
+
+  Xb_a = X * (beta .* s_a);
+  Xb_m = X * (beta .* s_m);
 }
 
 model {
@@ -119,5 +126,5 @@ model {
   beta ~ normal(0, sigmas);
 
   // Likelihood
-  y ~ normal(trend + X * beta, sigma_obs);
+  y ~ normal(trend .* (1 + Xb_m) + Xb_a, sigma_obs);
 }
diff --git a/python/stan/win/prophet.stan b/python/stan/win/prophet.stan
index 34e3163..a968ded 100644
--- a/python/stan/win/prophet.stan
+++ b/python/stan/win/prophet.stan
@@ -59,7 +59,7 @@ functions {
 
     gamma = logistic_gamma(k, m, delta, t_change, S);
     for (i in 1:T) {
-      Y[i] = cap[i] / (1 + exp(-(k + dot_product(A[i], delta)) * (t[i] - (m + dot_product(A[i], gamma)))))
+      Y[i] = cap[i] / (1 + exp(-(k + dot_product(A[i], delta)) * (t[i] - (m + dot_product(A[i], gamma)))));
     }
     return Y;
   }
@@ -81,7 +81,7 @@ functions {
     
     gamma = (-t_change .* delta);
     for (i in 1:T) {
-      Y[i] = (k + dot_product(A[i], delta)) * t[i] + (m + dot_product(A[i], gamma))
+      Y[i] = (k + dot_product(A[i], delta)) * t[i] + (m + dot_product(A[i], gamma));
     }
     return Y;
   }
@@ -99,6 +99,8 @@ data {
   vector[K] sigmas;     // Scale on seasonality prior
   real<lower=0> tau;    // Scale on changepoints prior
   int trend_indicator;  // 0 for linear, 1 for logistic
+  vector[K] s_a;        // Indicator of additive features
+  vector[K] s_m;        // Indicator of multiplicative features
 }
 
 transformed data {
@@ -117,6 +119,8 @@ parameters {
 transformed parameters {
   vector[T] trend;
   vector[T] Y;
+  vector[T] Xb_a;
+  vector[T] Xb_m;
 
   if (trend_indicator == 0) {
     trend = linear_trend(k, m, delta, t, A, t_change, S, T);
@@ -125,7 +129,9 @@ transformed parameters {
   }
 
   for (i in 1:T) {
-    Y[i] = trend[i] + dot_product(X[i], beta);
+    Xb_a[i] = dot_product(X[i], beta .* s_a);
+    Xb_m[i] = dot_product(X[i], beta .* s_m);
+    Y[i] = trend[i] * (1 + Xb_m[i]) + Xb_a[i];
   }
 }