We saw previously that when linear models are fitted to time series
data, the residuals can display autocorrelation.
When this occurs, it suggests that the model errors are not
independent.
In this lecture we shall examine:
- the consequences of failing to account for dependence between
errors;
- autoregressive models for a serially dependent sequence of
errors.
Consequences of Ignoring Autocorrelation in the Errors
Usually we fit a linear model to data by the method of least
squares.
Recall that under standard assumptions (A1 - A4), the least squares
estimators are unbiased for the regression parameters. In other words,
\(E[\hat \beta_i] = \beta_i\) for
i=1,…,p.
This unbiasedness property holds even if the error terms are
autocorrelated.
However, standard errors given by standard least squares theory (and
hence by standard R output) will usually be smaller than the true
standard deviations of the parameter estimates.
The result is that t-test statistics will be inflated, and
confidence intervals will be too narrow.
Accounting for Autocorrelation in the Errors
The previous discussion indicates that we need to account for
autocorrelation in the errors when it is present.
This will require two steps:
- Definition of an appropriate model for autocorrelated errors;
- Fitting the linear model with autocorrelated errors.
We will turn first to modelling for autocorrelated errors using
autoregressive models.
Autoregressive Models
In a time series, the errors \(\varepsilon_1, \varepsilon_2, \ldots,
\varepsilon_n\) can be thought of as a realization of a random
process (in time).
We can model serial correlation between these errors by making each
one explicitly dependent on the previous one: \[\varepsilon_t = \phi \varepsilon_{t-1} +
\eta_t\] where \(\eta_t~(t=1,2,\ldots,n)\) is white
noise, that is, a sequence of independent and identically
distributed (normal) random variables with zero mean.
- This equation is a first order autoregressive (or
AR1) model for \(\varepsilon_t\).
- Under this model, \(\mbox{Corr}(\varepsilon_t, \varepsilon_{t-1}) =
\phi\).
- More generally, \(\mbox{Corr}(\varepsilon_t, \varepsilon_{t-k}) =
\phi^k\)
ACF Plots for AR1 Processes
To illustrate an ACF plot for a function demonstrating
autocorrelation, we create two new variables from white noise
Z1 = rnorm(1000)
X1 = Z1[-1] + 0.85 * Z1[-1000]
X2 = Z1[-1] - 0.85 * Z1[-1000]
par(mfrow = c(2, 2))
plot(X1, type = "l")
acf(X1)
plot(X2, type = "l")
acf(X2)
The \(X1\) series has correlation of
0.467 at lag1, while the \(X2\) series
has autocorrelation of -0.448 at lag1
Autoregressive Models for Errors in Linear Modelling
If residuals in a linear model for time series data display
autocorrelation then we change the model to allow for an AR1 correlation
structure for the residuals in our model to account for this.
N.B. More complex models are also possible. For example:
- AR models of order q>1, where \(\varepsilon_t\) is modelled as explicitly
dependent on the previous q errors;
- moving average (MA) models and autoregressive moving average (ARMA)
models.
- For the purposes of this course we will focus on AR1 models for the
error process.
Application to the Tourism Data
Download Motel.csv
## Tourism <- read.csv(file = "Motel.csv", header = TRUE)
Tourism$Yr = Tourism$Year - 1979 + Tourism$Month/12
Tourism$Month <- factor(Tourism$Month)
Tourism.lm <- lm(RoomNights ~ Yr + Month, data = Tourism)
acf(resid(Tourism.lm))
There is a suggestion of exponential decay in autocorrelations. - AR1
model (with positive \(\phi\)) for
errors seems reasonable.
The model incorporating the AR1 error process, and written in
regression format, is as follows:
\[\begin{aligned}
Y_{t} &=& \beta_0 + \beta_1 \mbox{Yr} + \beta_2 z_{t,2} +
\beta_3 z_{t,3}
+ \ldots + \beta_{12} z_{t,12} + \varepsilon_t \\
\varepsilon_t &=& \phi \varepsilon_{t-1} +
\eta_{t}\end{aligned}\]
where:
- Yt is the tth observation of
RoomNights
;
- zt,j = 1 if observation t is taken in
month j;
- zt,j = 0 otherwise.
- \(\eta_1, \ldots, \eta_n\) are
independent and identically distributed normal random variables with
mean zero.
Model Fitting with Autoregressive Errors
To fit a model with AR1 errors, we must:
- Estimate the size of the autocorrelation parameter, \(\phi\).
- Account for the autocorrelation in errors when estimating the
regression coefficients.
The AR1 parameter, \(\phi\), can be
estimated by the method of (restricted) maximum likelihood.
The details are beyond the scope of this paper, but you will meet
maximum likelihood estimation if you take course 161.304.
Once \(\phi\) is estimated, we can
account for the AR1 nature of the errors by fitting the linear model
using generalized least squares.
Generalized Least Squares in R
The mathematical details of generalized least squares will not be
covered in this paper.
In R, linear models can be fitted using this technique using the
gls()
command found in the nlme
add-on
package. Call this package using library(nlme)
before using
the gls()
command.
- The
gls()
command uses the same model formula as
lm()
, but also accepts an additional argument to describe
the correlation structure for the errors.
Generalized Least Squares and the Tourism Data
Download Motel.csv
## Tourism <- read.csv(file = "Motel.csv", header = TRUE)
Tourism$Yr = Tourism$Year - 1979 + Tourism$Month/12
Tourism$Month <- factor(Tourism$Month)
Tourism.gls <- gls(RoomNights ~ Yr + Month, correlation = corAR1(), data = Tourism)
anova(Tourism.gls)
Denom. DF: 167
numDF F-value p-value
(Intercept) 1 45239.15 <.0001
Yr 1 1092.85 <.0001
Month 11 81.94 <.0001
Generalized least squares fit by REML
Model: RoomNights ~ Yr + Month
Data: Tourism
AIC BIC logLik
3722.941 3769.711 -1846.471
Correlation Structure: AR(1)
Formula: ~1
Parameter estimate(s):
Phi
0.392721
Coefficients:
Value Std.Error t-value p-value
(Intercept) 261682.96 4946.534 52.90229 0.0000
Yr 12747.00 382.435 33.33115 0.0000
Month2 -31194.09 4205.044 -7.41825 0.0000
Month3 23482.73 4959.531 4.73487 0.0000
Month4 -2170.47 5225.494 -0.41536 0.6784
Month5 -19419.12 5325.710 -3.64630 0.0004
Month6 -65827.11 5362.987 -12.27434 0.0000
Month7 -44522.34 5373.207 -8.28599 0.0000
Month8 -29768.74 5365.559 -5.54812 0.0000
Month9 -12039.27 5331.948 -2.25795 0.0252
Month10 26393.95 5238.620 5.03834 0.0000
Month11 16471.67 4988.022 3.30224 0.0012
Month12 -58631.94 4272.401 -13.72342 0.0000
Correlation:
(Intr) Yr Month2 Month3 Month4 Month5 Month6 Month7 Month8 Month9
Yr -0.640
Month2 -0.424 0.003
Month3 -0.498 -0.001 0.589
Month4 -0.522 -0.005 0.462 0.659
Month5 -0.529 -0.011 0.416 0.535 0.682
Month6 -0.529 -0.016 0.398 0.487 0.560 0.690
Month7 -0.526 -0.022 0.391 0.468 0.512 0.568 0.692
Month8 -0.521 -0.028 0.387 0.460 0.492 0.520 0.571 0.692
Month9 -0.514 -0.034 0.383 0.454 0.481 0.498 0.520 0.569 0.690
Month10 -0.501 -0.039 0.376 0.445 0.470 0.482 0.493 0.513 0.561 0.683
Month11 -0.473 -0.044 0.358 0.423 0.446 0.456 0.463 0.471 0.491 0.538
Month12 -0.400 -0.048 0.308 0.364 0.384 0.392 0.396 0.400 0.407 0.425
Mnth10 Mnth11
Yr
Month2
Month3
Month4
Month5
Month6
Month7
Month8
Month9
Month10
Month11 0.662
Month12 0.471 0.597
Standardized residuals:
Min Q1 Med Q3 Max
-2.5656600 -0.6285951 -0.1307054 0.5424996 3.2823396
Residual standard error: 14799.38
Degrees of freedom: 180 total; 167 residual
LS0tDQp0aXRsZTogIkxlY3R1cmUgMzU6IE1vZGVscyB3aXRoIEF1dG9yZWdyZXNzaXZlIEVycm9ycyINCnN1YnRpdGxlOiAxNjEuMjUxIFJlZ3Jlc3Npb24gTW9kZWxsaW5nDQphdXRob3I6ICJQcmVzZW50ZWQgYnkgSm9uYXRoYW4gTWFyc2hhbGwgPEouQy5tYXJzaGFsbEBtYXNzZXkuYWMubno+IiAgDQpkYXRlOiAiV2VlayAxMiBvZiBTZW1lc3RlciAyLCBgciBsdWJyaWRhdGU6OnllYXIobHVicmlkYXRlOjpub3coKSlgIg0Kb3V0cHV0Og0KICBodG1sX2RvY3VtZW50Og0KICAgIGNvZGVfZG93bmxvYWQ6IHRydWUNCiAgICB0aGVtZTogeWV0aQ0KICAgIGhpZ2hsaWdodDogdGFuZ28NCiAgaHRtbF9ub3RlYm9vazoNCiAgICBjb2RlX2Rvd25sb2FkOiB0cnVlDQogICAgdGhlbWU6IHlldGkNCiAgICBoaWdobGlnaHQ6IHRhbmdvDQogIGJlYW1lcl9wcmVzZW50YXRpb246IGRlZmF1bHQNCiAgd29yZF9kb2N1bWVudDogZGVmYXVsdA0KICBpb3NsaWRlc19wcmVzZW50YXRpb246IGRlZmF1bHQNCiAgc2xpZHlfcHJlc2VudGF0aW9uOiANCiAgICB0aGVtZTogeWV0aQ0KICAgIGhpZ2hsaWdodDogdGFuZ28NCiAgcGRmX2RvY3VtZW50OiBkZWZhdWx0DQotLS0NCg0KDQoNCg0KPCEtLS0gRGF0YSBpcyBvbg0KaHR0cHM6Ly9yLXJlc291cmNlcy5tYXNzZXkuYWMubnovZGF0YS8xNjEyNTEvDQotLS0+DQoNCmBgYHtyIHNldHVwLCBwdXJsPUZBTFNFLCBpbmNsdWRlPUZBTFNFfQ0KbGlicmFyeShrbml0cikNCm9wdHNfY2h1bmskc2V0KGRldj1jKCJwbmciLCAicGRmIikpDQpvcHRzX2NodW5rJHNldChmaWcuaGVpZ2h0PTYsIGZpZy53aWR0aD03LCBmaWcucGF0aD0iRmlndXJlcy8iLCBmaWcuYWx0PSJ1bmxhYmVsbGVkIikNCm9wdHNfY2h1bmskc2V0KGNvbW1lbnQ9IiIsIGZpZy5hbGlnbj0iY2VudGVyIiwgdGlkeT1UUlVFKQ0Kb3B0aW9ucyhrbml0ci5rYWJsZS5OQSA9ICcnKQ0KbGlicmFyeSh0aWR5dmVyc2UpDQpsaWJyYXJ5KGJyb29tKQ0KYGBgDQoNCg0KPCEtLS0gRG8gbm90IGVkaXQgYW55dGhpbmcgYWJvdmUgdGhpcyBsaW5lLiAtLS0+DQoNCg0KYGBge3IgZXh0cmFMaWIsIGluY2x1ZGU9RkFMU0V9IA0KbGlicmFyeShubG1lKSAgICAgICAgDQpgYGANCg0KDQpXZSBzYXcgcHJldmlvdXNseSB0aGF0IHdoZW4gbGluZWFyIG1vZGVscyBhcmUgZml0dGVkIHRvIHRpbWUgc2VyaWVzDQogICAgZGF0YSwgdGhlIHJlc2lkdWFscyBjYW4gZGlzcGxheSBhdXRvY29ycmVsYXRpb24uDQoNCldoZW4gdGhpcyBvY2N1cnMsIGl0IHN1Z2dlc3RzIHRoYXQgdGhlIG1vZGVsIGVycm9ycyBhcmUgbm90ICAgIGluZGVwZW5kZW50Lg0KDQpJbiB0aGlzIGxlY3R1cmUgd2Ugc2hhbGwgZXhhbWluZToNCg0KLSB0aGUgY29uc2VxdWVuY2VzIG9mIGZhaWxpbmcgdG8gYWNjb3VudCBmb3IgZGVwZW5kZW5jZSBiZXR3ZWVuDQogICAgICAgIGVycm9yczsNCi0gYXV0b3JlZ3Jlc3NpdmUgbW9kZWxzIGZvciBhIHNlcmlhbGx5IGRlcGVuZGVudCBzZXF1ZW5jZSBvZg0KICAgICAgICBlcnJvcnMuDQoNCiMjIENvbnNlcXVlbmNlcyBvZiBJZ25vcmluZyBBdXRvY29ycmVsYXRpb24gaW4gdGhlIEVycm9ycw0KDQpVc3VhbGx5IHdlIGZpdCBhIGxpbmVhciBtb2RlbCB0byBkYXRhIGJ5IHRoZSBtZXRob2Qgb2YgbGVhc3QNCiAgICBzcXVhcmVzLg0KDQpSZWNhbGwgdGhhdCB1bmRlciBzdGFuZGFyZCBhc3N1bXB0aW9ucyAoQTEgLSBBNCksIHRoZSBsZWFzdCBzcXVhcmVzIGVzdGltYXRvcnMNCiAgICBhcmUgdW5iaWFzZWQgZm9yIHRoZSByZWdyZXNzaW9uIHBhcmFtZXRlcnMuIEluIG90aGVyIHdvcmRzLCAkRVtcaGF0IFxiZXRhX2ldID0gXGJldGFfaSQgZm9yICppPTEsLi4uLHAqLg0KDQpUaGlzIHVuYmlhc2VkbmVzcyBwcm9wZXJ0eSBob2xkcyBldmVuIGlmIHRoZSBlcnJvciB0ZXJtcyBhcmUNCiAgICBhdXRvY29ycmVsYXRlZC4NCg0KSG93ZXZlciwgc3RhbmRhcmQgZXJyb3JzIGdpdmVuIGJ5IHN0YW5kYXJkIGxlYXN0IHNxdWFyZXMgdGhlb3J5IChhbmQgICAgIGhlbmNlIGJ5IHN0YW5kYXJkIFIgb3V0cHV0KSB3aWxsIHVzdWFsbHkgYmUgc21hbGxlciB0aGFuIHRoZSB0cnVlICAgICBzdGFuZGFyZCBkZXZpYXRpb25zIG9mIHRoZSBwYXJhbWV0ZXIgZXN0aW1hdGVzLg0KDQpUaGUgcmVzdWx0IGlzIHRoYXQgKnQqLXRlc3Qgc3RhdGlzdGljcyB3aWxsIGJlIGluZmxhdGVkLCBhbmQNCiAgICBjb25maWRlbmNlIGludGVydmFscyB3aWxsIGJlIHRvbyBuYXJyb3cuDQoNCiMjIEFjY291bnRpbmcgZm9yIEF1dG9jb3JyZWxhdGlvbiBpbiB0aGUgRXJyb3JzDQoNClRoZSBwcmV2aW91cyBkaXNjdXNzaW9uIGluZGljYXRlcyB0aGF0IHdlIG5lZWQgdG8gYWNjb3VudCBmb3INCiAgICBhdXRvY29ycmVsYXRpb24gaW4gdGhlIGVycm9ycyB3aGVuIGl0IGlzIHByZXNlbnQuDQoNClRoaXMgd2lsbCByZXF1aXJlIHR3byBzdGVwczoNCiAgICANCi0gRGVmaW5pdGlvbiBvZiBhbiBhcHByb3ByaWF0ZSBtb2RlbCBmb3IgYXV0b2NvcnJlbGF0ZWQgZXJyb3JzOw0KLSBGaXR0aW5nIHRoZSBsaW5lYXIgbW9kZWwgd2l0aCBhdXRvY29ycmVsYXRlZCBlcnJvcnMuDQoNCldlIHdpbGwgdHVybiBmaXJzdCB0byBtb2RlbGxpbmcgZm9yIGF1dG9jb3JyZWxhdGVkIGVycm9ycyB1c2luZyAgICAgYXV0b3JlZ3Jlc3NpdmUgbW9kZWxzLg0KDQojIyBBdXRvcmVncmVzc2l2ZSBNb2RlbHMNCg0KSW4gYSB0aW1lIHNlcmllcywgdGhlIGVycm9ycyAkXHZhcmVwc2lsb25fMSwgXHZhcmVwc2lsb25fMiwgXGxkb3RzLCBcdmFyZXBzaWxvbl9uJA0KICAgIGNhbiBiZSB0aG91Z2h0IG9mIGFzIGEgcmVhbGl6YXRpb24gb2YgYSByYW5kb20gcHJvY2VzcyAoaW4gdGltZSkuDQoNCldlIGNhbiBtb2RlbCBzZXJpYWwgY29ycmVsYXRpb24gYmV0d2VlbiB0aGVzZSBlcnJvcnMgYnkgbWFraW5nIGVhY2gNCiAgICBvbmUgZXhwbGljaXRseSBkZXBlbmRlbnQgb24gdGhlIHByZXZpb3VzIG9uZToNCiAgICAkJFx2YXJlcHNpbG9uX3QgPSBccGhpIFx2YXJlcHNpbG9uX3t0LTF9ICsgXGV0YV90JCQgd2hlcmUgJFxldGFfdH4odD0xLDIsXGxkb3RzLG4pJCBpcyAqKndoaXRlICBub2lzZSoqLCB0aGF0IGlzLCBhIHNlcXVlbmNlIG9mIGluZGVwZW5kZW50IGFuZCBpZGVudGljYWxseSBkaXN0cmlidXRlZCAgKG5vcm1hbCkgcmFuZG9tIHZhcmlhYmxlcyB3aXRoIHplcm8gbWVhbi4NCg0KLSBUaGlzIGVxdWF0aW9uICBpcyBhIGZpcnN0IG9yZGVyIGF1dG9yZWdyZXNzaXZlIChvciAgICAgKipBUjEqKikgbW9kZWwgZm9yICRcdmFyZXBzaWxvbl90JC4NCi0gVW5kZXIgdGhpcyBtb2RlbCwgJFxtYm94e0NvcnJ9KFx2YXJlcHNpbG9uX3QsIFx2YXJlcHNpbG9uX3t0LTF9KSA9IFxwaGkkLg0KLSBNb3JlIGdlbmVyYWxseSwgJFxtYm94e0NvcnJ9KFx2YXJlcHNpbG9uX3QsIFx2YXJlcHNpbG9uX3t0LWt9KSA9IFxwaGleayQNCg0KIyMgQUNGIFBsb3RzIGZvciBBUjEgUHJvY2Vzc2VzDQoNClRvIGlsbHVzdHJhdGUgYW4gQUNGIHBsb3QgZm9yIGEgZnVuY3Rpb24gZGVtb25zdHJhdGluZyBhdXRvY29ycmVsYXRpb24sIHdlIGNyZWF0ZSB0d28gbmV3IHZhcmlhYmxlcyBmcm9tIHdoaXRlIG5vaXNlIA0KDQpgYGB7ciB3aGl0ZU5vaXNlLCBmaWcud2lkdGg9OCwgZmlnLmhlaWdodD04fQ0KWjEgPSBybm9ybSgxMDAwKQ0KWDEgPSBaMVstMV0gKyAwLjg1KloxWy0xMDAwXQ0KWDIgPSBaMVstMV0gLSAwLjg1KloxWy0xMDAwXQ0KcGFyKG1mcm93PWMoMiwyKSkNCnBsb3QoWDEsIHR5cGU9ImwiKQ0KYWNmKFgxKQ0KcGxvdChYMiwgdHlwZT0ibCIpDQphY2YoWDIpDQpgYGANCg0KVGhlICRYMSQgc2VyaWVzIGhhcyBjb3JyZWxhdGlvbiBvZiBgciByb3VuZChjb3IoWDFbLTFdLFgxWy05OTldKSwzKWAgYXQgbGFnMSwgd2hpbGUgdGhlICRYMiQgc2VyaWVzIGhhcyBhdXRvY29ycmVsYXRpb24gb2YgYHIgcm91bmQoY29yKFgyWy0xXSxYMlstOTk5XSksMylgIGF0IGxhZzENCg0KDQojIyBBdXRvcmVncmVzc2l2ZSBNb2RlbHMgZm9yIEVycm9ycyBpbiBMaW5lYXIgTW9kZWxsaW5nDQoNCklmIHJlc2lkdWFscyBpbiBhIGxpbmVhciBtb2RlbCBmb3IgdGltZSBzZXJpZXMgZGF0YSBkaXNwbGF5DQogICAgYXV0b2NvcnJlbGF0aW9uIHRoZW4gd2UgY2hhbmdlIHRoZSBtb2RlbCB0byBhbGxvdyBmb3IgYW4gQVIxIGNvcnJlbGF0aW9uIHN0cnVjdHVyZSBmb3IgdGhlIHJlc2lkdWFscyBpbiBvdXIgbW9kZWwgdG8gYWNjb3VudCBmb3IgdGhpcy4NCg0KTi5CLiBNb3JlIGNvbXBsZXggbW9kZWxzIGFyZSBhbHNvIHBvc3NpYmxlLiBGb3IgZXhhbXBsZToNCiAgICANCi0gQVIgbW9kZWxzIG9mIG9yZGVyICpxPjEqLCB3aGVyZSAkXHZhcmVwc2lsb25fdCQgaXMgbW9kZWxsZWQgYXMNCiAgICAgICAgZXhwbGljaXRseSBkZXBlbmRlbnQgb24gdGhlIHByZXZpb3VzICpxKiBlcnJvcnM7DQotIG1vdmluZyBhdmVyYWdlIChNQSkgbW9kZWxzIGFuZCBhdXRvcmVncmVzc2l2ZSBtb3ZpbmcgYXZlcmFnZQ0KICAgICAgICAoQVJNQSkgbW9kZWxzLg0KLSBGb3IgdGhlIHB1cnBvc2VzIG9mIHRoaXMgY291cnNlIHdlIHdpbGwgZm9jdXMgb24gQVIxIG1vZGVscyBmb3IgdGhlDQogICAgZXJyb3IgcHJvY2Vzcy4NCg0KIyMgQXBwbGljYXRpb24gdG8gdGhlIFRvdXJpc20gRGF0YQ0KYHIgeGZ1bjo6ZW1iZWRfZmlsZSgiLi4vLi4vZGF0YS9Nb3RlbC5jc3YiKWANCg0KDQpgYGB7ciBnZXRUb3VyaXNtRGF0YSwgZWNobz0tMSwgZXZhbD0tMn0NClRvdXJpc20gPC0gcmVhZC5jc3YoZmlsZT0iLi4vLi4vZGF0YS9Nb3RlbC5jc3YiLCBoZWFkZXI9VFJVRSkNClRvdXJpc20gPC0gcmVhZC5jc3YoZmlsZT0iTW90ZWwuY3N2IiwgaGVhZGVyPVRSVUUpDQpUb3VyaXNtJFlyID0gVG91cmlzbSRZZWFyIC0xOTc5ICsgVG91cmlzbSRNb250aC8xMg0KVG91cmlzbSRNb250aCA8LSBmYWN0b3IoVG91cmlzbSRNb250aCkNCmBgYA0KDQpgYGB7ciBUb3VyaXNtLmxtUmVzaWRBQ0Z9DQpUb3VyaXNtLmxtIDwtIGxtKFJvb21OaWdodHN+WXIgKyBNb250aCwgZGF0YT1Ub3VyaXNtKQ0KYWNmKHJlc2lkKFRvdXJpc20ubG0pKQ0KYGBgDQoNClRoZXJlIGlzIGEgc3VnZ2VzdGlvbiBvZiBleHBvbmVudGlhbCBkZWNheSBpbiBhdXRvY29ycmVsYXRpb25zLg0KLSBBUjEgbW9kZWwgKHdpdGggcG9zaXRpdmUgJFxwaGkkKSBmb3IgZXJyb3JzIHNlZW1zIHJlYXNvbmFibGUuDQoNCg0KVGhlIG1vZGVsIGluY29ycG9yYXRpbmcgdGhlIEFSMSBlcnJvciBwcm9jZXNzLCBhbmQgd3JpdHRlbiBpbiByZWdyZXNzaW9uIGZvcm1hdCwNCmlzIGFzIGZvbGxvd3M6IA0KDQokJFxiZWdpbnthbGlnbmVkfQ0KWV97dH0gJj0mIFxiZXRhXzAgKyBcYmV0YV8xIFxtYm94e1lyfSArIFxiZXRhXzIgel97dCwyfSArIFxiZXRhXzMgel97dCwzfSANCisgXGxkb3RzICsgXGJldGFfezEyfSB6X3t0LDEyfSArIFx2YXJlcHNpbG9uX3QgXFwNClx2YXJlcHNpbG9uX3QgJj0mIFxwaGkgXHZhcmVwc2lsb25fe3QtMX0gKyBcZXRhX3t0fVxlbmR7YWxpZ25lZH0kJCANCg0Kd2hlcmU6DQoNCi0gKll+dH4qIGlzIHRoZSAqdCpedGheIG9ic2VydmF0aW9uIG9mIGBSb29tTmlnaHRzYDsNCi0gKnp+dCxqfiA9IDEqIGlmIG9ic2VydmF0aW9uICp0KiBpcyB0YWtlbiBpbiBtb250aCAqaio7DQotICp6fnQsan4gPSAwKiBvdGhlcndpc2UuDQotICRcZXRhXzEsIFxsZG90cywgXGV0YV9uJCBhcmUgaW5kZXBlbmRlbnQgYW5kIGlkZW50aWNhbGx5DQogICAgZGlzdHJpYnV0ZWQgbm9ybWFsIHJhbmRvbSB2YXJpYWJsZXMgd2l0aCBtZWFuIHplcm8uDQoNCiMjIE1vZGVsIEZpdHRpbmcgd2l0aCBBdXRvcmVncmVzc2l2ZSBFcnJvcnMNCg0KVG8gZml0IGEgbW9kZWwgd2l0aCBBUjEgZXJyb3JzLCB3ZSBtdXN0Og0KDQotIEVzdGltYXRlIHRoZSBzaXplIG9mIHRoZSBhdXRvY29ycmVsYXRpb24gcGFyYW1ldGVyLCAkXHBoaSQuDQotIEFjY291bnQgZm9yIHRoZSBhdXRvY29ycmVsYXRpb24gaW4gZXJyb3JzIHdoZW4gZXN0aW1hdGluZyB0aGUNCiAgICAgICAgcmVncmVzc2lvbiBjb2VmZmljaWVudHMuDQoNClRoZSBBUjEgcGFyYW1ldGVyLCAkXHBoaSQsIGNhbiBiZSBlc3RpbWF0ZWQgYnkgdGhlIG1ldGhvZCBvZg0KICAgIChyZXN0cmljdGVkKSBtYXhpbXVtIGxpa2VsaWhvb2QuDQoNClRoZSBkZXRhaWxzIGFyZSBiZXlvbmQgdGhlIHNjb3BlIG9mIHRoaXMgcGFwZXIsIGJ1dCB5b3Ugd2lsbCBtZWV0ICAgICBtYXhpbXVtIGxpa2VsaWhvb2QgZXN0aW1hdGlvbiBpZiB5b3UgdGFrZSBjb3Vyc2UgMTYxLjMwNC4NCg0KT25jZSAkXHBoaSQgaXMgZXN0aW1hdGVkLCB3ZSBjYW4gYWNjb3VudCBmb3IgdGhlIEFSMSBuYXR1cmUgb2YgdGhlDQogICAgZXJyb3JzIGJ5IGZpdHRpbmcgdGhlIGxpbmVhciBtb2RlbCB1c2luZyAqKmdlbmVyYWxpemVkIGxlYXN0DQogICAgc3F1YXJlcyoqLg0KDQojIyBHZW5lcmFsaXplZCBMZWFzdCBTcXVhcmVzIGluIFINCg0KVGhlIG1hdGhlbWF0aWNhbCBkZXRhaWxzIG9mIGdlbmVyYWxpemVkIGxlYXN0IHNxdWFyZXMgd2lsbCBub3QgYmUgICAgIGNvdmVyZWQgaW4gdGhpcyBwYXBlci4NCg0KSW4gUiwgbGluZWFyIG1vZGVscyBjYW4gYmUgZml0dGVkIHVzaW5nIHRoaXMgdGVjaG5pcXVlIHVzaW5nIHRoZQ0KICAgIGBnbHMoKWAgY29tbWFuZCBmb3VuZCBpbiB0aGUgYG5sbWVgIGFkZC1vbiBwYWNrYWdlLiBDYWxsIHRoaXMgcGFja2FnZSB1c2luZyBgbGlicmFyeShubG1lKWAgYmVmb3JlIHVzaW5nIHRoZSBgZ2xzKClgIGNvbW1hbmQuDQoNCg0KICAgIA0KDQpgYGB7ciBnZXRMaWJ9IA0KbGlicmFyeShubG1lKSAgICAgICAgDQpgYGANCg0KLSBUaGUgYGdscygpYCBjb21tYW5kIHVzZXMgdGhlIHNhbWUgbW9kZWwgZm9ybXVsYSBhcyBgbG0oKWAsIGJ1dCBhbHNvIGFjY2VwdHMgYW4gYWRkaXRpb25hbCBhcmd1bWVudCB0byBkZXNjcmliZSAgICB0aGUgY29ycmVsYXRpb24gc3RydWN0dXJlIGZvciB0aGUgZXJyb3JzLg0KDQojIyBHZW5lcmFsaXplZCBMZWFzdCBTcXVhcmVzIGFuZCB0aGUgVG91cmlzbSBEYXRhDQoNCmByIHhmdW46OmVtYmVkX2ZpbGUoIi4uLy4uL2RhdGEvTW90ZWwuY3N2IilgDQoNCmBgYHtyIFRvdXJpc20uZ2xzLCBlY2hvPS0xLCBldmFsPS0yfQ0KVG91cmlzbSA8LSByZWFkLmNzdihmaWxlPSIuLi8uLi9kYXRhL21vdGVsLmNzdiIsIGhlYWRlcj1UUlVFKQ0KVG91cmlzbSA8LSByZWFkLmNzdihmaWxlPSJNb3RlbC5jc3YiLCBoZWFkZXI9VFJVRSkNClRvdXJpc20kWXIgPSBUb3VyaXNtJFllYXIgLTE5NzkgKyBUb3VyaXNtJE1vbnRoLzEyDQpUb3VyaXNtJE1vbnRoIDwtIGZhY3RvcihUb3VyaXNtJE1vbnRoKQ0KVG91cmlzbS5nbHMgPC0gZ2xzKFJvb21OaWdodHN+WXIgKyBNb250aCwgY29ycmVsYXRpb249Y29yQVIxKCksIGRhdGE9VG91cmlzbSkgDQphbm92YShUb3VyaXNtLmdscykNCmBgYA0KDQoNCmBgYHtyIHN1bW1hcnlUb3VyaXNtLmdsc30NCnN1bW1hcnkoVG91cmlzbS5nbHMpDQpgYGANCg0KDQojIyMgQ29tbWVudHMgb24gdGhlIEFuYWx5c2lzDQoNClRoZSB1c2Ugb2YgYW4gQVIxIGVycm9yIHByb2Nlc3MgaXMgb2J0YWluZWQgYnkgc2V0dGluZyAgDQogICAgYGNvcnJlbGF0aW9uPWNvckFSMSgpYCAgDQogICAgaW4gdGhlIGBnbHMoKWAgY29tbWFuZC4NCg0KRXJyb3IgdmFyaWFuY2UgYW5kIGNvcnJlbGF0aW9uIGVzdGltYXRlZCBieSBSRU1MIChyZXN0cmljdGVkDQogICAgbWF4aW11bSBsaWtlbGlob29kKS4NCg0KRXN0aW1hdGVkIHZhbHVlIG9mIEFSMSBjb3JyZWxhdGlvbiBwYXJhbWV0ZXIgaXMNCiAgICAkXGhhdCBccGhpID0gMC4zOTMkLg0KDQpCb3RoIGBZcmAgYW5kIGBNb250aGAgYXJlIGFzc29jaWF0ZWQgd2l0aA0KICAgIFJvb21OaWdodHMgKCpQKi12YWx1ZSBsZXNzIHRoYW4gJDAuMDAwMSQgaW4gYm90aA0KICAgIGNhc2VzKS4NCg0KQ29tcGFyaXNvbiBvZiB0aGUgb3V0cHV0IGZyb20gYHN1bW1hcnkoVG91cmlzbS5nbHMpYCB3aXRoDQogICAgdGhlIGNvcnJlc3BvbmRpbmcgb3V0cHV0IGZyb20gdGhlIGVhcmxpZXIgYW5hbHlzaXMgKGFzc3VtaW5nDQogICAgaW5kZXBlbmRlbnQgZXJyb3JzKSBzaG93cyB0aGF0IHRoZSBzdGFuZGFyZCBlcnJvcnMgYXJlIHVzdWFsbHkNCiAgICBsYXJnZXIgaGVyZSwgcGFydGljdWxhcmx5IGZvciB0aGUgaW50ZXJjZXB0IGFuZCB0aGUgY29lZmZpY2llbnQgb2YgICAgYFlyYC4NCg0KVGhlIG9ubHkgcmVtYWluaW5nIHF1ZXN0aW9uIGlzLCB3YXMgdGhlIHVzZSBvZiBgZ2xzKClgIHN1Y2Nlc3NmdWw/DQoNCmBgYHtyIGdldFJlc2lkcywgaW5jbHVkZT1UUlVFfQ0KYWNmKHJlc2lkdWFscyhUb3VyaXNtLmdscywgdHlwZT0ibiIpKQ0KDQpgYGANCg0KV2hhdCBtaWdodCBiZSB0aGUgbmV4dCBzdGVwPyAgSXQgbG9va3MgbGlrZSB0aGVyZSBpcyBhIHNpZ25pZmljYW50IDEyLW1vbnRoIGNvcnJlbGF0aW9uLCB3aGljaCB3b3VsZCBhbHNvIG1ha2Ugc2Vuc2UuIA0KDQpXZSB3aWxsIGhhdmUgdG8gbGVhdmUgdGhlIHRhc2sgb2YgYWRqdXN0aW5nIGZvciBtb3JlIGNvbXBsaWNhdGVkIGNvcnJlbGF0aW9uIHN0cnVjdHVyZXMgZm9yIGFub3RoZXIgY291cnNlLg0KDQoNCg==
Comments on the Analysis
The use of an AR1 error process is obtained by setting
correlation=corAR1()
in the
gls()
command.Error variance and correlation estimated by REML (restricted maximum likelihood).
Estimated value of AR1 correlation parameter is \(\hat \phi = 0.393\).
Both
Yr
andMonth
are associated with RoomNights (P-value less than \(0.0001\) in both cases).Comparison of the output from
summary(Tourism.gls)
with the corresponding output from the earlier analysis (assuming independent errors) shows that the standard errors are usually larger here, particularly for the intercept and the coefficient ofYr
.The only remaining question is, was the use of
gls()
successful?What might be the next step? It looks like there is a significant 12-month correlation, which would also make sense.
We will have to leave the task of adjusting for more complicated correlation structures for another course.