# How to implement model in R?

I would like your help to implement this model in R

or more explicity

where

• yt = monthly mean values
• μi = mean value in month i, i = 1 . . . 12 .
• I1;t = Indicator series for month i of the year, i.e., 1 if the month corresponds to month i of the year, and 0 otherwise.
• bi = Trend in month i of the year.
• Rt = t/12
• Z1;t= Regressor 1, with c1 the associated coefficient.
• Z2;t= Regressor 2, with c2 the associated coefficient.
• Nt = Residual noise series, modeled as an autoregressive AR(1) series

so i would like to get 12 coefficients for monthly and trend component and 1 for regressors 1 and 2.

Here is one way to fit the model that you describe.

# sample series
x <- AirPassengers
# to illustrate a more general case,
# take a subsample that does not start in the first season
# and ends in the last season
x <- window(x, start=c(1949,2), end=c(1959,4))


Indicator variables for the seasonal intercepts can be created in several ways, for example:

# monthly intercepts
S <- frequency(x)
monthly.means <- do.call("rbind",
replicate(ceiling(length(x)/S), diag(S), simplify = FALSE))


Some arrangements for non-square time series (this has no effect if the series starts in the first season and ends in the last season):

monthly.means <- ts(monthly.means, frequency = S, start = c(start(x)[1], 1))
monthly.means <- window(monthly.means, start = start(x), end = end(x))


Regardless of the month the first observation belongs to, the first column is related to January, the second to February and so on.

# column names
if (S == 12) {
colnames(monthly.means) <- month.abb
} else
colnames(monthly.means) <- paste("season", 1L:S)


First rows of the monthly.means object:

monthly.means[1:15,]
#       Jan Feb Mar Apr May Jun Jul Aug Sep Oct Nov Dec
#  [1,]   0   1   0   0   0   0   0   0   0   0   0   0
#  [2,]   0   0   1   0   0   0   0   0   0   0   0   0
#  [3,]   0   0   0   1   0   0   0   0   0   0   0   0
#  [4,]   0   0   0   0   1   0   0   0   0   0   0   0
#  [5,]   0   0   0   0   0   1   0   0   0   0   0   0
#  [6,]   0   0   0   0   0   0   1   0   0   0   0   0
#  [7,]   0   0   0   0   0   0   0   1   0   0   0   0
#  [8,]   0   0   0   0   0   0   0   0   1   0   0   0
#  [9,]   0   0   0   0   0   0   0   0   0   1   0   0
# [10,]   0   0   0   0   0   0   0   0   0   0   1   0
# [11,]   0   0   0   0   0   0   0   0   0   0   0   1
# [12,]   1   0   0   0   0   0   0   0   0   0   0   0
# [13,]   0   1   0   0   0   0   0   0   0   0   0   0
# [14,]   0   0   1   0   0   0   0   0   0   0   0   0
# [15,]   0   0   0   1   0   0   0   0   0   0   0   0


For the monthly trends we can reuse monthly.means':

monthly.trends <- monthly.means * seq_along(x) / S
round(monthly.trends[1:15,], 2)
#       Jan  Feb  Mar  Apr  May  Jun Jul  Aug  Sep  Oct  Nov  Dec
#  [1,]   0 0.08 0.00 0.00 0.00 0.00 0.0 0.00 0.00 0.00 0.00 0.00
#  [2,]   0 0.00 0.17 0.00 0.00 0.00 0.0 0.00 0.00 0.00 0.00 0.00
#  [3,]   0 0.00 0.00 0.25 0.00 0.00 0.0 0.00 0.00 0.00 0.00 0.00
#  [4,]   0 0.00 0.00 0.00 0.33 0.00 0.0 0.00 0.00 0.00 0.00 0.00
#  [5,]   0 0.00 0.00 0.00 0.00 0.42 0.0 0.00 0.00 0.00 0.00 0.00
#  [6,]   0 0.00 0.00 0.00 0.00 0.00 0.5 0.00 0.00 0.00 0.00 0.00
#  [7,]   0 0.00 0.00 0.00 0.00 0.00 0.0 0.58 0.00 0.00 0.00 0.00
#  [8,]   0 0.00 0.00 0.00 0.00 0.00 0.0 0.00 0.67 0.00 0.00 0.00
#  [9,]   0 0.00 0.00 0.00 0.00 0.00 0.0 0.00 0.00 0.75 0.00 0.00
# [10,]   0 0.00 0.00 0.00 0.00 0.00 0.0 0.00 0.00 0.00 0.83 0.00
# [11,]   0 0.00 0.00 0.00 0.00 0.00 0.0 0.00 0.00 0.00 0.00 0.92
# [12,]   1 0.00 0.00 0.00 0.00 0.00 0.0 0.00 0.00 0.00 0.00 0.00
# [13,]   0 1.08 0.00 0.00 0.00 0.00 0.0 0.00 0.00 0.00 0.00 0.00
# [14,]   0 0.00 1.17 0.00 0.00 0.00 0.0 0.00 0.00 0.00 0.00 0.00
# [15,]   0 0.00 0.00 1.25 0.00 0.00 0.0 0.00 0.00 0.00 0.00 0.00


Some arbitrary regressors:

set.seed(123)
xreg1 <- runif(length(x), 100, 200)
xreg2 <- rnorm(length(x), mean = mean(x))


Edited

As I did in my first edit, the model could now be fitted using the function lm.

lm(x ~ 0 + monthly.means + monthly.trends + xreg1 + xreg2)


But I overlooked the AR(1) structure for the error term that is mentioned in the question. As @RobHyndman says in the comment below, we can use the function arima specifying an AR(1) process for the error term, order = c(1,0,0). An intercept should not be included in order to avoid exact multicollinearity with the dummies for the seasonal means. (The example is arbitrary, the output is printed just to show the coefficients that are estimated.)

xreg <- cbind(monthly.means, monthly.trends, xreg1, xreg2)
fit <- arima(x, order = c(1,0,0), xreg = xreg, include.mean = FALSE)
round(coef(fit), 2)
#                ar1  monthly.means.Jan  monthly.means.Feb  monthly.means.Mar
#               0.84              69.87              82.45              93.10
#  monthly.means.Apr  monthly.means.May  monthly.means.Jun  monthly.means.Jul
#              85.17              76.39              77.48              83.31
#  monthly.means.Aug  monthly.means.Sep  monthly.means.Oct  monthly.means.Nov
#              81.01              81.23              68.06              57.34
#  monthly.means.Dec monthly.trends.Jan monthly.trends.Feb monthly.trends.Mar
#              73.20              27.07              23.36              27.84
# monthly.trends.Apr monthly.trends.May monthly.trends.Jun monthly.trends.Jul
#              27.56              29.34              35.98              40.47
# monthly.trends.Aug monthly.trends.Sep monthly.trends.Oct monthly.trends.Nov
#              40.30              32.07              27.85              24.02
# monthly.trends.Dec              xreg1              xreg2
#              25.60               0.00               0.08


We can do some simple check. In a model containing only the seasonal intercepts, the estimates of the coefficients match the values of the sample means.

fit2 <- arima(x, order = c(0,0,0), xreg = xreg[,1:12], include.mean = FALSE)
coef.seasonal.intercepts <- coef(fit2)
sample.means <- lapply(split(x, cycle(x)), mean)
all.equal(coef.seasonal.intercepts,
as.vector(unlist(sample.means)), check.attributes = FALSE)
# TRUE


Edited (Answer to a comment posted by the OP.)

If the disturbance term follows an AR(1) process, the Ordinary Least Squares estimator is unbiased but not efficient, i.e. on average it gives the true value but the standard errors of parameters estimates are higher than in the classical setting of independent errors (in other words, estimates are not efficient).

As you say, extending the model with the AR(1) error term via arima will not change much the the estimates (standard OLS is still unbiased), but their standard errors will be smaller due to the gain in efficiency.

When the disturbance term is autocorrelated, omitting the AR term will lead to larger standard errors for the estimates, which implies a larger denominator in the t-statistics and therefore larger t-statistics. Hence, the tests for the null of non-significant regressors will be biased towards non-rejection. Using arima to specify the AR error term will protect against this issue.

The code below is a small exercise to check these ideas: 10,000 series are generated from a model with intercept and AR(1) errors. If xregcoef is set to a value other than zero, then a external variable is added to the data. A model that includes the intercept and the external regressor is fitted by means of arima with order=c(1,0,0) (stored in fit1) and by means of lm (stored in fit2).

set.seed(123)
xreg <- runif(200, 2, 6)
xregcoef <- 0
res <- matrix(nrow = 10000, ncol = 6)
colnames(res) <- c("coef 1", "s.e. 1", "t-stat. 1", "coef 2", "s.e. 2", "t-stat. 2")
for (i in seq.int(nrow(res)))
{
x <- 2 + arima.sim(n=200, model=list(ar=0.7)) + xregcoef * xreg
fit1 <- arima(x, order=c(1,0,0), xreg=xreg, include.mean=TRUE)
res[i,1] <- coef(fit1)["xreg"]
res[i,2] <- sqrt(fit1\$var.coef["xreg","xreg"])
res[i,3] <- res[i,1]/res[i,2]
fit2 <- summary(lm(x ~ 1 + xreg))
res[i,4:6] <- coef(fit2)["xreg", c("Estimate", "Std. Error", "t value")]
}


Setting for example xregcoef=3, we can see that estimates from both regressions are very close to the true value (unbiased estimates) but the standard errors are slightly higher for lm, where when the AR structure is omitted.

# results for xregcoef=3
t(apply(res, 2, summary))
#               Min.  1st Qu.   Median     Mean  3rd Qu.     Max.
# coef 1     2.81600  2.96300  3.00100  3.00000  3.03800  3.20100
# s.e. 1     0.04313  0.05286  0.05491  0.05492  0.05693  0.06677
# coef 2     2.65800  2.93300  3.00100  2.99900  3.06400  3.37000
# s.e. 2     0.05928  0.08393  0.08879  0.08905  0.09400  0.12280


Let's now consider xregcoef=0, that is, the regressor is not part of the data generating process but we fit a model including this regressor. On average the coefficient is estimated as zero in both cases, however, as the standard errors are slightly higher in the second case, the null of the t-test is rejected slightly more often than it should, given a 5% significance level chosen below.

# results for xregcoef=0
t(apply(res, 2, summary))
#               Min.  1st Qu.    Median       Mean 3rd Qu.    Max.
# coef 1    -0.18440 -0.03698 0.0007062  0.0003874 0.03763 0.20080
# s.e. 1     0.04313  0.05286 0.0549100  0.0549200 0.05693 0.06677
# t-stat. 1 -3.68100 -0.67410 0.0129200  0.0072410 0.68500 3.53400
# coef 2    -0.34200 -0.06694 0.0012510 -0.0006191 0.06422 0.36970
# s.e. 2     0.05928  0.08393 0.0887900  0.0890500 0.09400 0.12280
# t-stat. 2 -3.96600 -0.75560 0.0138800 -0.0070140 0.72210 4.44500
#
# rejections of the null xregcoef=0 at the 5% significance level
sum(abs(res[,3]) > 1.96) / nrow(res)
# [1] 0.0516
sum(abs(res[,6]) > 1.96) / nrow(res)
# [1] 0.0749
`
• Nice work, but you miss the AR(1) error. Use cbind() to collect the regressors into X, then arima() with xreg=X, order=c(1,0,0) and include.mean=FALSE. – Rob Hyndman Nov 23 '14 at 21:38
• @RobHyndman I fixed it, many thanks. – javlacalle Nov 23 '14 at 23:31
• Thank you very much for your attention and your time, I aprrecite it very much. – Giovana Miranda Nov 24 '14 at 17:59
• Please correct me if i am wrong, but fitting the arima model instead of the multilinear, wont affect much the coefficients values, but will change their statistical signicance and their standard errors? – Giovana Miranda Nov 24 '14 at 23:19
• @GiovanaMiranda I have added some content in response to your comment. – javlacalle Nov 25 '14 at 12:04