Is there a way to allow seasonality in regression coefficients? Say I have a time series, Gt, and a covariate Bt. I want to find the relationship between them by the ARMA model: 
Gt = Zt + β0 + β1Bt
where the residual Zt follows some ARMA process. 
The problem is: I know for sure that β0 and β1 varies with the time of the year. Yet I do not want to fit a separate model to each month because that introduces discontinuity into my time series, which means I cannot calculate the autocorrelation function of the final residuals. 
So, is there a time series model (or family of models, I wonder) that allows the correlation coefficients of its covariates to change seasonally? 
========================
Edit: Thank you for those who replied here. I decided to just use seasonal dummies, but got busy so failed to reply in time. 
 A: Edit
(The same idea was proposed by Stephan Kolassa a few minutes before I posted my answer. The answer below can still give you some relevant details.)
You could use seasonal dummies. For simplicity I illustrate this for a quarterly time series. Seasonal dummies are indicator variables for each season. The $i$-th seasonal dummy takes on the value 1 for those observations related to season $i$ and 0 otherwise. For a quarterly series the seasonal dummies, $SD$, are defined as follows:
\begin{eqnarray}
SD = \left[
\begin{array}{cccc}
1 & 0 & 0 & 0 \\
0 & 1 & 0 & 0 \\
0 & 0 & 1 & 0 \\
0 & 0 & 0 & 1 \\
1 & 0 & 0 & 0 \\
\vdots & \vdots & \vdots & \vdots \\
1 & 0 & 0 & 0 \\
0 & 1 & 0 & 0 \\
0 & 0 & 1 & 0 \\
0 & 0 & 0 & 1 \\
\end{array}
\right] \quad
SDB = \left[
\begin{array}{cccc}
B_{1} & 0 & 0 & 0 \\
0 & B_{2} & 0 & 0 \\
0 & 0 & B_{3} & 0 \\
0 & 0 & 0 & B_{4} \\
B_{5} & 0 & 0 & 0 \\
\vdots & \vdots & \vdots & \vdots \\
B_{n-3} & 0 & 0 & 0 \\
0 & B_{n-2} & 0 & 0 \\
0 & 0 & B_{n-1} & 0 \\
0 & 0 & 0 & B_{n} \\
\end{array}
\right]
\end{eqnarray}
You can multiply each column in $SD$ by your explanatory variable $B_t$ and get the matrix $SDB$ defined above.
Then, you can specify your model as follows:
$$
G_t = Z_t + \beta_{0,s} SD_t + \beta_{1,s} SDB_t \,,
$$
where the index $s$ indicates the season. Observe that we now have four coefficients (12 in your monthly series) $\beta_{1,s}$, one for each column in $SDB$. 
The same for the intercept $\beta_0$ except that we must remove one column in $SD$ in order to avoid perfect collinearity. In a monthly series you would include for example the first 11 seasonal intercepts in $SD$.
Fitting the model for example by maximum likelihood will give you one coefficient estimate for each season. You could also test whether $\beta_{0,s}$ are the same for all $s$ or similarly if $\beta_{1,s}$ are constant across seasons.
A: Certainly there is. Simply include monthly dummies in an interaction with $B_t$. Let $M_{tm}$ denote a dummy that is 1 if time $t$ corresponds to month $m$ and 0 otherwise. Then fit the following regression with ARMA errors:
$$ G_t = \beta M_{t\cdot} + \gamma B_tM_{t\cdot} + Z_t $$
where $Z_t$ is ARMA(p,q) and $\beta$ and $\gamma$ are parameter vectors of length 12.
You can do the actual fitting using R with the nlme package, using the gls() function and specifying a corARMA() correlation structure.
A: If you don't want to discretise the seasonal effect, you could assume that the regression coefficients vary in a cyclic manner as a function of the time of year, i.e. $\beta_0(t) = w_0 + w_1\sin nt + w_2\cos nt$ and $\beta_1(t) = w_3 + w_4\sin nt + w_5\cos nt$, then if you substitute these into your linear model, you should get something of the form
$G_t = Z_t + w_o +  w_1\sin nt + w_2\cos nt + w_3B_t + w_4B_t\sin nt + w_5B_t\cos nt$
You could fit this model by using OLS regression (or whatever method you are already using) with the additional covariates $\sin nt$, $\cos nt$, $B_t\sin nt$ and $B_t\cos nt$, where $n$ is whatever constant you need to represent a year ($2\pi/365$ for a daily time-series).
This wouldn't introduce any discontinuities in the model as the seasonality in the regression coefficients are smooth functions of time.  I suspect if you added sine and cosine components representing harmonics of the annual cycle you could model deviations from simple sinusoidal variation in the regression coefficients (Fourier series type approach).
Caveat: Been a long day, so I may have made a stupid error somewhere.
A: Fit the mean and the harmonics of the seasonal cycle to the time series of x and y. These provide the intercept terms. Then, subtract them from x and y to create anomalies. Use these anomalies x' and y' to compute seasonally varying regression slope coefficients: Fit the array product between the x' and y' with the mean and leading harmonics to the seasonal cycle. Do the same for the variance of the x'. Then divide the seasonal cycle fit to the covariance by the seasonal cycle fit to the variance to provide continuously evolving slope coefficients. For details, see
http://onlinelibrary.wiley.com/doi/10.1002/qj.3054/full
