Time Series Analysis - Spring 2018,
Cyrus Safaie
March 30, 2018
1 Week1 Overview
1.1 Agenda
Introduction to Time Series
White noise and random walk
Correlation, independence, orthogonality and auto-correlation
Stationarity and non-stationarity
1.2 What is forecastable?
how well we understand the factors that contribute to it;
how much data are available;
whether the forecasts can affect the thing we are trying to forecast.
1.3 What to forecast?
granularity: every product line, or for groups of products?
granularity: every sales outlet, or for outlets grouped by region, or only for total sales?
weekly data, monthly data or annual data?
1.4 Forecasting Types:
- Cross-sectional forecasting: Cross-sectional models are used when the variable to be forecast exhibits a relationship with one or more other predictor variables. Models in this class include regression models, additive models, and some kinds of neural networks
Some people use the term “predict” for cross-sectional data and “forecast” for time series data
- Time series forecasting:Time series data are useful when you are forecasting something that is changing over time. Anything that is observed sequentially over time is a time series. Stock price, whether. When forecasting time series data, the aim is to estimate how the sequence of observations will continue into the future
Time series forecasting uses only information on the variable to be forecast, and makes no attempt to discover the factors which affect its behavior. Therefore it will extrapolate trend and seasonal patterns, but it ignores all other information
you can also set up time series with predictors. This is where a lot of new techniques are being added
\(ED=f(current temperature, strength of economy, population,time of day, day of week, error).\)
OR
\(ED_{t+1}=f(ED_{t},ED_{t−1},ED_{t−2},ED_{t−3},…,error)\)
1.5 Autocorrelation
correlation measures the extent of a linear relationship between two variables, auto-correlation measures the linear relationship between lagged values of a time series.
There are several auto-correlation coefficients, depending on the lag length
\(r_{k}= \dfrac{\sum_{t=k+1}^{T}(y_{t}-\bar{y})(y_{t-k}-\bar{y})}{\sum_{t=1}^{T} (y_{t}-\bar{y})}\)
T is the length of the time series.
1.6 ACF
If you plot the auto-correlation coefficients you get the ACF plot
Example from textbook
suppressWarnings(library(forecast))
suppressWarnings(library(fpp))## Loading required package: fma## Loading required package: expsmooth## Loading required package: lmtest## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric## Loading required package: tseriesbeer2 <- window(ausbeer, start=1992, end=2006-.1)
lag.plot(beer2, lags=9, do.lines=FALSE)
head(beer2)## Qtr1 Qtr2 Qtr3 Qtr4
## 1992 443 410 420 532
## 1993 433 421Acf(beer2)
1.7 White noise:
Time series that show no auto-correlation are called “white noise” (auto-correlation close to zero). There are none zero due to randomness.
For a white noise series, we expect 95% of the spikes in the ACF to lie within \(\dfrac{±2} {\sqrt{T}}\) where T is the length of the time series. That’s the blue line on ACF plot
set.seed(60661)
x <- ts(rnorm(150))
plot(x, main="White noise")
acf(x)
1.8 Random Walk
Mathematically represented as \(Y_{t}=Y_{t-1} + w_{t}\)
Autocovariance–measures the linear dependency between two points in the same TS observed at different times.
1.9 Stationarity
A stationary time series is one whose properties do not depend on the time at which the series is observed.
In mathematics and statistics, a stationary process (a.k.a. a strict(ly) stationary process or strong(ly) stationary process) is a stochastic process whose unconditional joint probability distribution does not change when shifted in time. Consequently, parameters such as mean and variance, if they are present, also do not change over time.
if \(y_{t}\) is a stationary time series, then for all \(s\), the distribution of \((y_{t},…,y_{t+s})\) does not depend on \(t\)
Basic idea is that the laws of probability that govern the process behavior do not change over time.
STRICTLY stationary –when the probability behaviors of every collection of values of a TS is identical to that of the time shifted T {s+ℎ}for all h
WEAKLY stationary -is a finite variance process such that:
- Mean is constant and does not depend on
- Autocovariance function depends on depend on t and s through their absolute delta
- If TS is Gaussian then weak stationarity=== strict stationarity
So time series with trends, or with seasonality, are not stationary — the trend and seasonality will affect the value of the time series at different times.
White noise series is stationary — it does not matter when you observe it, it should look much the same at any period of time.
a stationary time series will have no predictable patterns in the long-term.
For a stationary time series, the ACF will drop to zero relatively quickly, while the ACF of non-stationary data decreases slowly.
for non-stationary data, the value of r1r1 is often large and positive.
adf test and kpss tests for stationarity
1.9.1 Differencing
Differencing can help stabilize the mean of a time series by removing changes in the level of a time series, and so eliminating trend and seasonality.
2 Week2 Overview
2.1 Agenda
Smoothing Time Series (Moving Average / Exponential)
Holt Winters
Regression analysis
Univariate and multivariate regression modeling – model assumptions & multicollinearity
2.2 Simple and multiple Regression (again…)
The “error” term does not imply a mistake, but a deviation from the underlying straight line model. It captures anything that may affect yiyi other than xixi. We assume that these errors:
have mean zero; otherwise the forecasts will be systematically biased.
are not auto-correlated; otherwise the forecasts will be inefficient as there is more information to be exploited in the data.
are unrelated to the predictor variable; otherwise there would be more information that should be included in the systematic part of the model.
It is also useful to have the errors normally distributed with constant variance in order to produce prediction intervals and to perform statistical inference. While these additional conditions make the calculations simpler, they are not necessary for forecasting.
The forecast values of y obtained from the observed x values are called “fitted values”.
Rememebr Residuals is same as error
2.2.1 Regression and correlation
\(r\) measures the strength and the direction (positive or negative) of the linear relationship between the two variables. The stronger the linear relationship, the closer the observed data points will cluster around a straight line.
\(\hat{B_{1}} = r * S_{y}/S_{x}\)
where s is standard deviation
2.2.2 residual review:
A non-random pattern may indicate that a non-linear relationship may be required
or some heteroscedasticity is present (i.e., the residuals show non-constant variance)
or there is some left over serial correlation (only when the data are time series).
2.2.3 Non-linear functional forms
\(log (y_{i})= \beta_{0} + \beta_{1}log (x_{i})+ \epsilon_{i}\)
In this model, the slope β1β1 can be interpeted as an elasticity: β1β1 is the average percentage change in yy resulting from a 1%1% change in xx.
Example from text book:
par(mfrow=c(1,2))
fit2 <- lm(log(Carbon) ~ log(City), data=fuel)
plot(jitter(Carbon) ~ jitter(City), xlab="City (mpg)",
ylab="Carbon footprint (tonnes per year)", data=fuel)
lines(1:50, exp(fit2$coef[1]+fit2$coef[2]*log(1:50)))
plot(log(jitter(Carbon)) ~ log(jitter(City)),
xlab="log City mpg", ylab="log carbon footprint", data=fuel)
abline(fit2)
2.2.4 Regression with time series data
fit.ex3 <- tslm(consumption ~ income, data=usconsumption)
plot(usconsumption, ylab="% change in consumption and income",
plot.type="single", col=1:2, xlab="Year")
legend("topright", legend=c("Consumption","Income"),
lty=1, col=c(1,2), cex=.9)
plot(consumption ~ income, data=usconsumption,
ylab="% change in consumption", xlab="% change in income")
abline(fit.ex3)
summary(fit.ex3)##
## Call:
## tslm(formula = consumption ~ income, data = usconsumption)
##
## Residuals:
## Min 1Q Median 3Q Max
## -2.3681 -0.3237 0.0266 0.3436 1.5581
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 0.52062 0.06231 8.356 2.79e-14 ***
## income 0.31866 0.05226 6.098 7.61e-09 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 0.6274 on 162 degrees of freedom
## Multiple R-squared: 0.1867, Adjusted R-squared: 0.1817
## F-statistic: 37.18 on 1 and 162 DF, p-value: 7.614e-09Challenges:
future values of the predictor variable (Income in this case) are needed to be input into the estimated model, but these are not known in advance. often you have to fortecast those separately
when fitting a regression model to time series data, it is very common to find autocorrelation in the residuals.It is not wrong becasue it’s not biased
Spurious regression More often than not, time series data are “non-stationary”; that is, the values of the time series do not fluctuate around a constant mean or with a constant variance. High \(R_{2}\) and high residual autocorrelation can be signs of spurious regression. often when there are highly correlated data that are totally unrelated such as “#of air passengers”, and “Rice production”. Great example in text book. This often given good short term forecast but fails long term
beer2 <- window(ausbeer,start=1992,end=2006-.1)
fit <- tslm(beer2 ~ trend + season)
summary(fit)##
## Call:
## tslm(formula = beer2 ~ trend + season)
##
## Residuals:
## Min 1Q Median 3Q Max
## -44.024 -8.390 0.249 8.619 23.320
##
## Coefficients:
## Estimate Std. Error t value Pr(>|t|)
## (Intercept) 441.8141 4.5338 97.449 < 2e-16 ***
## trend -0.3820 0.1078 -3.544 0.000854 ***
## season2 -34.0466 4.9174 -6.924 7.18e-09 ***
## season3 -18.0931 4.9209 -3.677 0.000568 ***
## season4 76.0746 4.9268 15.441 < 2e-16 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
## Residual standard error: 13.01 on 51 degrees of freedom
## Multiple R-squared: 0.921, Adjusted R-squared: 0.9149
## F-statistic: 148.7 on 4 and 51 DF, p-value: < 2.2e-162.2.5 Multicollinearity
A closely related issue is multicollinearity which occurs when similar information is provided by two or more of the predictor variables in a multiple regression. It can occur in a number of ways.
Two predictors are highly correlated with each other
A linear combination of predictors is highly correlated with another linear combination of predictors
In most statistical softwares, if you are not interested in the specific contributions of each predictor, and if the future values of your predictor variables are within their historical ranges, there is nothing to worry about — multicollinearity is not a problem.
2.3 Time series decomp
- Trend
A trend exists when there is a long-term increase or decrease in the data. It does not have to be linear. Sometimes we will refer to a trend “changing direction” when it might go from an increasing trend to a decreasing trend.
- Seasonal
A seasonal pattern exists when a series is influenced by seasonal factors (e.g., the quarter of the year, the month, or day of the week). Seasonality is always of a fixed and known period.
- Cyclic
A cyclic pattern exists when data exhibit rises and falls that are not of fixed period.
2.3.1 additive model
most appropriate if the magnitude of the seasonal fluctuations or the variation around the trend-cycle does not vary with the level of the time series.
\(y_{t}=S_{t}+T_{t}+E_{t}\)
e.g. toy sales increase by $1 million every Dec.
where \(y_{t}\) is the data at period t, St is the seasonal component at period t, \(T_{t}\) is the trend-cycle component at period t and Et is the remainder (or irregular or error) component at period t.
2.3.2 multiplicative model (common with econometrics models)
When the variation in the seasonal pattern, or the variation around the trend-cycle, appears to be proportional to the level of the time series, then a multiplicative model is more appropriate.
\(y_{t}=S_{t}T_{t}E_{t}\)
OR
\(log(y_{t})=log(S_{t})+log(T_{t})+log(E_{t})\)
e.g. toy sales increase by 42% every Dec.
example from book
fit <- stl(elecequip, s.window=5)
plot(elecequip, col="gray",
main="Electrical equipment manufacturing",
ylab="New orders index", xlab="")
lines(fit$time.series[,2],col="red",ylab="Trend")
plot(fit)
Forecasts of STL objects are obtained by applying a non-seasonal forecasting method to the seasonally adjusted data and re-seasonalizing using the last year of the seasonal component.
- For an additive model, the seasonally adjusted data are given by yt−St, and for multiplicative data, the seasonally adjusted values are obtained using yt/St
2.4 moving average
smoother:
\(\hat{T}_{t}= 1/m \sum_{j=-k}^{K}(y_{t+j})\)
where \(m=2k+1\). That is, the estimate of the trend-cycle at time t is obtained by averaging values of the time series within k periods of t.
2.5 Exponential smoothing:
2.5.1 simple exponential smoothing
This method is suitable for forecasting data with no trend or seasonal pattern.
wighted is what we talk about in the lecture that more recent observation have more weight
2.5.2 Holt’s linear trend method
Extended simple exponential smoothing to allow forecasting of data with a trend.
2.5.3 Holt-Winters seasonal method
Exponential smoothing with trend and seasonality The Holt-Winters seasonal method comprises the forecast equation and three smoothing equations
3 Week3 Overview
3.1 Agenda
Data transformations
Box-Jenkins ARMA models
Box-Jenkins ARIMA models
Stationarity and invertibility
Model specification
3.2 ARIMA Model
Exponential smoothing models are based on a description of trend and seasonality in the data, ARIMA models aim to describe the autocorrelations in the data.
3.2.1 Stationarity
A stationary time series is one whose properties do not depend on the time at which the series is observed. More precisely, if \(y_{t}\) is a stationary time series, then for all \(s\), the distribution of \((y_{t},…,y_{t+s})\) does not depend on \(t\).
Time series with trends, or with seasonality, are not stationary
White noise series is stationary — it does not matter when you observe it
A time series with cyclic behaviour (but not trend or seasonality) is also stationary. That is because the cycles are not of fixed length,
In general, a stationary time series will have no predictable patterns in the long-term.
- ACF plots are a good way identify non stationary time series
How to determine stationarity?
- For a stationary time series, the ACF will drop to zero relatively quickly, while the ACF of non-stationary data decreases slowly.
** One way to determine more objectively if differencing is required is to use a unit root test. ADF (\(H_{0}\) is non-stationarity), KPSS, etc.
3.2.2 Differencing
Differencing is a way to make a time series stationary by computing the differences between consecutive observations.
Transformations such as logarithms can help to stabilize the variance of a time series.
Differencing can help stabilize the mean of a time series by removing changes in the level of a time series, and so eliminating trend and seasonality.
3.2.3 Random walk models:
When the differenced series is white noise, the model can be written as:
\(y_{t}−y_{t-1}=e_{t}\)
where \(e\) is white noise
- Widely used for non-stationary data especially in finance and econ
Random walks typically have:
- long periods of apparent trends up or down
- sudden and unpredictable changes in direction.
The forecasts from a random walk model are equal to the last observation, as future movements are unpredictable
- sometime the differncing is vs last year same time (seasonal naive) or vs more than one previous period
3.2.4 Autoregressive models (AR(p))
In an autoregression model, we forecast the variable of interest using a linear combination of past values of the variable.
\(y_{t}=c+ϕ1y_{t-1}1+ϕ2y_{t-2}+⋯+ϕp y_{t-p}+e_{t}\),
When ϕ1=0, yt is equivalent to white noise.
When ϕ1=1 and c=0, yt is equivalent to a random walk.
When ϕ1=1 and c≠0, yt is equivalent to a random walk with drift
3.2.5 Moving average models (MA(q))
Rather than use past values of the forecast variable in a regression, a moving average model uses past forecast errors in a regression-like model. yt can be thought of as a weighted moving average of the past few forecast errors
\(y_{t}=c+e_{t} + ϕ1 e_{t-1} + ϕ 2e_{t-2}+⋯+ϕ pe_{t-q}\)
- Note: A moving average model is used for forecasting future values while moving average smoothing is used for estimating the trend-cycle of past values.
3.2.6 Non-seasonal ARIMA models
If we combine differencing with autoregression and a moving average model, we obtain a non-seasonal ARIMA model.The “predictors” on the right hand side include both lagged values of ytyt and lagged errors. We call this an ARIMA(p,d,qp,d,q) model, where
p = order of the autoregressive part;
d = degree of first differencing involved for \(y_{t}\) ;
q = order of the moving average part.
finding p, d, and q are not trivial but with auto.arima it makes it pretty easy, but with automation you need to be careful with d and c to ensure things are not breaking and you are also searching enough space to find AIC optimum
White noise ARIMA(0,0,0)
Random walk ARIMA(0,1,0) with no constant
Random walk with drift ARIMA(0,1,0) with a constant
Autoregression ARIMA(p,0,0)
Moving average ARIMA(0,0,q)
3.2.6.1 PACF
It is sometimes possible to use the ACF plot, and the closely related PACF plot, to determine appropriate values for p and q as well.
- partial autocorrelations. These measure the {relationship} between yt and yt−k after removing the effects of other time lags
par(mfrow=c(1,2))
Acf(usconsumption[,1],main="")
Pacf(usconsumption[,1],main="")
If the data are from an ARIMA(p,d,0p,d,0) or ARIMA(0,d,q0,d,q) model, then the ACF and PACF plots can be helpful in determining the value of p or q. If both p and q are positive, then the plots do not help in finding suitable values of p and q.
3.2.7 Choosing your own model
Plot the data. Identify any unusual observations.
If necessary, transform the data (using a Box-Cox transformation) to stabilize the variance.
If the data are non-stationary: take first differences of the data until the data are stationary.
Examine the ACF/PACF: Is an AR(pp) or MA(qq) model appropriate?
Try your chosen model(s), and use the AICc to search for a better model.
Check the residuals from your chosen model by plotting the ACF of the residuals, and doing a portmanteau test of the residuals. If they do not look like white noise, try a modified model.
Once the residuals look like white noise, calculate forecasts.
auto.arima only take care of step 3 to 5
4 Week4 Overview
4.1 Agenda
Akaike and Bayesian Information Criterion
Model estimation
Model diagnostics –residual analysis
Ljung Box test
Trend and seasonal decompositions
Seasonal ARIMA models
Final Project: Choose a time-series data applicationand develop a model solution
We will use course material for this week, most of this we have already discussed
4.1.1 Akaike and Bayesian Information Criterion
Hyndman recommend AICc for model selection
4.1.2 Modeling Procedure
Plot the data. Identify any unusual observations.
If necessary, transform the data (using a Box-Cox transformation) to stabilize the variance.
If the data are non-stationary: take first differences of the data until the data are stationary.library(fpp) elecequip
Examine the ACF/PACF: Is an AR(pp) or MA(qq) model appropriate?
Try your chosen model(s), and use the AICc to search for a better model.
Check the residuals from your chosen model by plotting the ACF of the residuals, and doing a portmanteau test of the residuals. If they do not look like white noise, try a modified model.
Once the residuals look like white noise, calculate forecasts.
auto.arima only take care of step 3 to 5
Here is a good example from the book we will review together:
library(fpp)
elecequip## Jan Feb Mar Apr May Jun Jul Aug Sep Oct
## 1996 79.43 75.86 86.40 72.67 74.93 83.88 79.88 62.47 85.50 83.19
## 1997 78.72 77.49 89.94 81.35 78.76 89.59 83.75 69.87 91.18 89.52
## 1998 81.97 85.26 93.09 81.19 85.74 91.24 83.56 66.45 93.45 86.03
## 1999 81.68 81.68 91.35 79.55 87.08 96.71 98.10 79.22 103.68 101.00
## 2000 95.42 98.49 116.37 101.09 104.20 114.79 107.75 96.23 123.65 116.24
## 2001 100.69 102.99 119.21 92.56 98.86 111.26 96.25 79.81 102.18 96.28
## 2002 89.66 89.23 104.36 87.17 89.43 102.25 88.26 75.73 99.60 96.57
## 2003 89.45 86.87 98.94 85.62 85.31 101.22 91.93 77.01 104.50 99.83
## 2004 89.93 92.73 105.22 91.56 92.60 104.46 96.28 79.61 105.55 99.15
## 2005 91.73 90.45 105.56 92.15 91.23 108.95 99.33 83.30 110.85 104.99
## 2006 99.09 99.73 116.05 103.51 102.99 119.45 107.98 90.50 121.85 117.12
## 2007 103.92 103.97 125.63 104.69 108.36 123.09 108.88 93.98 121.94 116.79
## 2008 109.35 105.64 121.30 108.62 103.13 117.84 103.62 89.22 109.41 103.93
## 2009 77.33 75.01 86.31 74.09 74.09 85.58 79.84 65.24 87.92 84.45
## 2010 79.16 78.40 94.32 84.45 84.92 103.18 89.42 77.66 95.68 94.03
## 2011 91.47 87.66 103.33 88.56 92.32 102.21 92.80 76.44 94.00 91.67
## Nov Dec
## 1996 84.29 89.79
## 1997 91.12 92.97
## 1998 86.91 93.42
## 1999 99.52 111.94
## 2000 117.00 128.75
## 2001 101.38 109.97
## 2002 96.22 101.12
## 2003 101.10 109.16
## 2004 99.81 113.72
## 2005 107.10 114.38
## 2006 113.66 120.35
## 2007 115.78 127.28
## 2008 100.07 101.15
## 2009 87.93 102.42
## 2010 100.99 101.26
## 2011 91.98plot(elecequip)
You can decompose a model as below to trend, seasonality, and stochastic parts
library(forecast)
eeadj2 <- stl(elecequip, s.window="periodic")
head(eeadj2)## $time.series
## seasonal trend remainder
## Jan 1996 -5.563267 80.49243 4.500836205
## Feb 1996 -6.088288 80.40279 1.545497118
## Mar 1996 7.977940 80.31315 -1.891091696
## Apr 1996 -6.347420 80.27340 -1.255982031
## May 1996 -4.818401 80.23365 -0.485249582
## Jun 1996 7.751458 80.22800 -4.099458302
## Jul 1996 -1.542437 80.22235 1.200086960
## Aug 1996 -16.876314 80.26219 -0.915875748
## Sep 1996 7.357938 80.30203 -2.159966935
## Oct 1996 2.931033 80.64139 -0.382421462
## Nov 1996 3.767879 80.98075 -0.458626817
## Dec 1996 11.449882 81.47659 -3.136471262
## Jan 1997 -5.563267 81.97243 2.310835875
## Feb 1997 -6.088288 82.48176 1.096531613
## Mar 1997 7.977940 82.99108 -1.029022376
## Apr 1997 -6.347420 83.48437 4.213049758
## May 1997 -4.818401 83.97766 -0.399255324
## Jun 1997 7.751458 84.36204 -2.523493620
## Jul 1997 -1.542437 84.74642 0.546022067
## Aug 1997 -16.876314 85.06688 1.679436299
## Sep 1997 7.357938 85.38734 -1.565277949
## Oct 1997 2.931033 85.69779 0.891179567
## Nov 1997 3.767879 86.00824 1.343886255
## Dec 1997 11.449882 86.19200 -4.671879363
## Jan 1998 -5.563267 86.37576 1.157506600
## Feb 1998 -6.088288 86.36255 4.985733976
## Mar 1998 7.977940 86.34935 -1.237288374
## Apr 1998 -6.347420 86.18547 1.351947663
## May 1998 -4.818401 86.02159 4.536806484
## Jun 1998 7.751458 85.80234 -2.313795990
## Jul 1998 -1.542437 85.58308 -0.480644482
## Aug 1998 -16.876314 85.36424 -2.037930105
## Sep 1998 7.357938 85.14541 0.946655792
## Oct 1998 2.931033 85.15061 -2.051642276
## Nov 1998 3.767879 85.15581 -2.013691171
## Dec 1998 11.449882 85.65811 -3.687990372
## Jan 1999 -5.563267 86.16040 1.082862010
## Feb 1999 -6.088288 87.16132 0.606966720
## Mar 1999 7.977940 88.16224 -4.790178296
## Apr 1999 -6.347420 89.37903 -3.481612189
## May 1999 -4.818401 90.59582 1.302576703
## Jun 1999 7.751458 91.94683 -2.988289328
## Jul 1999 -1.542437 93.29784 6.344598624
## Aug 1999 -16.876314 94.83100 1.265316097
## Sep 1999 7.357938 96.36416 -0.042094909
## Oct 1999 2.931033 97.96311 0.105861974
## Nov 1999 3.767879 99.56205 -3.809931971
## Dec 1999 11.449882 100.95387 -0.463752308
## Jan 2000 -5.563267 102.34569 -1.362421063
## Feb 2000 -6.088288 103.72084 0.857449840
## Mar 2000 7.977940 105.09599 3.296071016
## Apr 2000 -6.347420 106.52199 0.915432026
## May 2000 -4.818401 107.94799 1.070415820
## Jun 2000 7.751458 109.03172 -1.993174501
## Jul 2000 -1.542437 110.11545 -0.823010839
## Aug 2000 -16.876314 110.55891 2.547399435
## Sep 2000 7.357938 111.00238 5.289681229
## Oct 2000 2.931033 110.80663 2.502339414
## Nov 2000 3.767879 110.61087 2.621246771
## Dec 2000 11.449882 109.85816 7.441961963
## Jan 2001 -5.563267 109.10544 -2.852171264
## Feb 2001 -6.088288 107.74246 1.335831881
## Mar 2001 7.977940 106.37948 4.852585299
## Apr 2001 -6.347420 104.75838 -5.850961876
## May 2001 -4.818401 103.13729 0.541113734
## Jun 2001 7.751458 101.72695 1.781593868
## Jul 2001 -1.542437 100.31661 -2.524172016
## Aug 2001 -16.876314 99.25565 -2.569333652
## Sep 2001 7.357938 98.19469 -3.372623767
## Oct 2001 2.931033 97.41498 -4.066015638
## Nov 2001 3.767879 96.63528 0.976841663
## Dec 2001 11.449882 96.07237 2.447747812
## Jan 2002 -5.563267 95.50946 -0.286194456
## Feb 2002 -6.088288 95.11914 0.199152661
## Mar 2002 7.977940 94.72881 1.653250050
## Apr 2002 -6.347420 94.38498 -0.867560919
## May 2002 -4.818401 94.04115 0.207250897
## Jun 2002 7.751458 93.68108 0.817459362
## Jul 2002 -1.542437 93.32102 -3.518578192
## Aug 2002 -16.876314 93.04606 -0.439746619
## Sep 2002 7.357938 92.77111 -0.529043525
## Oct 2002 2.931033 92.57692 1.062050534
## Nov 2002 3.767879 92.38273 0.069393765
## Dec 2002 11.449882 92.33785 -2.667735075
## Jan 2003 -5.563267 92.29298 2.720287668
## Feb 2003 -6.088288 92.45669 0.501597219
## Mar 2003 7.977940 92.62040 -1.658342957
## Apr 2003 -6.347420 92.95235 -0.984930425
## May 2003 -4.818401 93.28430 -3.155895109
## Jun 2003 7.751458 93.72038 -0.251832828
## Jul 2003 -1.542437 94.15645 -0.684016566
## Aug 2003 -16.876314 94.64761 -0.761295578
## Sep 2003 7.357938 95.13877 2.003296930
## Oct 2003 2.931033 95.65165 1.247321701
## Nov 2003 3.767879 96.16453 1.167595645
## Dec 2003 11.449882 96.53822 1.171898584
## Jan 2004 -5.563267 96.91191 -1.418646896
## Feb 2004 -6.088288 97.07947 1.738819391
## Mar 2004 7.977940 97.24702 -0.004964048
## Apr 2004 -6.347420 97.29194 0.615483928
## May 2004 -4.818401 97.33685 0.081554690
## Jun 2004 7.751458 97.39972 -0.691178685
## Jul 2004 -1.542437 97.46260 0.359841923
## Aug 2004 -16.876314 97.50618 -1.019861883
## Sep 2004 7.357938 97.54976 0.642305831
## Oct 2004 2.931033 97.58748 -1.368515877
## Nov 2004 3.767879 97.62521 -1.583088414
## Dec 2004 11.449882 97.78707 4.483049285
## Jan 2005 -5.563267 97.94893 -0.655661434
## Feb 2005 -6.088288 98.26889 -1.730598621
## Mar 2005 7.977940 98.58885 -1.006785534
## Apr 2005 -6.347420 99.01145 -0.514028034
## May 2005 -4.818401 99.43405 -3.385647748
## Jun 2005 7.751458 99.94137 1.257172761
## Jul 2005 -1.542437 100.44869 0.423747252
## Aug 2005 -16.876314 101.16118 -0.984870672
## Sep 2005 7.357938 101.87368 1.618382925
## Oct 2005 2.931033 102.74467 -0.685697944
## Nov 2005 3.767879 103.61565 -0.283529641
## Dec 2005 11.449882 104.47443 -1.544311271
## Jan 2006 -5.563267 105.33321 -0.679941319
## Feb 2006 -6.088288 106.15247 -0.334180043
## Mar 2006 7.977940 106.97173 1.100331507
## Apr 2006 -6.347420 107.75508 2.102339581
## May 2006 -4.818401 108.53843 -0.730029560
## Jun 2006 7.751458 109.11715 2.581397265
## Jul 2006 -1.542437 109.69586 -0.173421928
## Aug 2006 -16.876314 110.10651 -2.730194049
## Sep 2006 7.357938 110.51716 3.974905351
## Oct 2006 2.931033 110.86354 3.325424999
## Nov 2006 3.767879 111.20993 -1.317806181
## Dec 2006 11.449882 111.47877 -2.578650106
## Jan 2007 -5.563267 111.74761 -2.264342450
## Feb 2007 -6.088288 111.91067 -1.852377772
## Mar 2007 7.977940 112.07372 5.578337180
## Apr 2007 -6.347420 112.24718 -1.209762541
## May 2007 -4.818401 112.42064 0.757760523
## Jun 2007 7.751458 112.69219 2.646349621
## Jul 2007 -1.542437 112.96374 -2.541307300
## Aug 2007 -16.876314 113.10571 -2.249399511
## Sep 2007 7.357938 113.24768 1.334379798
## Oct 2007 2.931033 113.18380 0.675163765
## Nov 2007 3.767879 113.11992 -1.107803095
## Dec 2007 11.449882 112.83168 2.998434357
## Jan 2008 -5.563267 112.54344 2.369823392
## Feb 2008 -6.088288 111.90263 -0.174339028
## Mar 2008 7.977940 111.26181 2.060248826
## Apr 2008 -6.347420 110.07213 4.895286277
## May 2008 -4.818401 108.88245 -0.934053487
## Jun 2008 7.751458 106.93074 3.157800638
## Jul 2008 -1.542437 104.97903 0.183408746
## Aug 2008 -16.876314 102.34177 3.754541933
## Sep 2008 7.357938 99.70452 2.347546640
## Oct 2008 2.931033 96.84375 4.155215883
## Nov 2008 3.767879 93.98299 2.319134298
## Dec 2008 11.449882 91.37215 -1.672034872
## Jan 2009 -5.563267 88.76132 -5.868052460
## Feb 2009 -6.088288 86.66833 -5.570043985
## Mar 2009 7.977940 84.57535 -6.243285237
## Apr 2009 -6.347420 83.31327 -2.875850398
## May 2009 -4.818401 82.05119 -3.142792774
## Jun 2009 7.751458 81.82624 -3.997698949
## Jul 2009 -1.542437 81.60129 -0.218851143
## Aug 2009 -16.876314 82.08196 0.034354607
## Sep 2009 7.357938 82.56263 -2.000568122
## Oct 2009 2.931033 83.39340 -1.874428133
## Nov 2009 3.767879 84.22416 -0.062038973
## Dec 2009 11.449882 85.25203 5.718091756
## Jan 2010 -5.563267 86.27989 -1.556625934
## Feb 2010 -6.088288 87.23455 -2.746259355
## Mar 2010 7.977940 88.18920 -1.847142502
## Apr 2010 -6.347420 88.97408 1.823341669
## May 2010 -4.818401 89.75895 -0.020551375
## Jun 2010 7.751458 90.42993 4.998614043
## Jul 2010 -1.542437 91.10090 -0.138466558
## Aug 2010 -16.876314 91.74133 2.794988946
## Sep 2010 7.357938 92.38175 -4.059684031
## Oct 2010 2.931033 92.85066 -1.751689222
## Nov 2010 3.767879 93.31957 3.902554760
## Dec 2010 11.449882 93.60016 -3.790040101
## Jan 2011 -5.563267 93.88075 3.152516621
## Feb 2011 -6.088288 93.78838 -0.040089325
## Mar 2011 7.977940 93.69601 1.656055003
## Apr 2011 -6.347420 93.29526 1.612159589
## May 2011 -4.818401 92.89451 4.243886960
## Jun 2011 7.751458 92.47931 1.979229392
## Jul 2011 -1.542437 92.06411 2.278325807
## Aug 2011 -16.876314 91.59660 1.719716179
## Sep 2011 7.357938 91.12908 -4.487021928
## Oct 2011 2.931033 90.58267 -1.843702891
## Nov 2011 3.767879 90.03626 -1.824134682
##
## $weights
## [1] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [36] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [71] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [106] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [141] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
## [176] 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1 1
##
## $call
## stl(x = elecequip, s.window = "periodic")
##
## $win
## s t l
## 1911 19 13
##
## $deg
## s t l
## 0 1 1
##
## $jump
## s t l
## 192 2 2plot(eeadj2)
Now we adjust for seasonality (basically remove the seasonality) for this example for
eeadj <- seasadj(stl(elecequip, s.window="periodic"))
eeadj## Jan Feb Mar Apr May Jun Jul
## 1996 84.99327 81.94829 78.42206 79.01742 79.74840 76.12854 81.42244
## 1997 84.28327 83.57829 81.96206 87.69742 83.57840 81.83854 85.29244
## 1998 87.53327 91.34829 85.11206 87.53742 90.55840 83.48854 85.10244
## 1999 87.24327 87.76829 83.37206 85.89742 91.89840 88.95854 99.64244
## 2000 100.98327 104.57829 108.39206 107.43742 109.01840 107.03854 109.29244
## 2001 106.25327 109.07829 111.23206 98.90742 103.67840 103.50854 97.79244
## 2002 95.22327 95.31829 96.38206 93.51742 94.24840 94.49854 89.80244
## 2003 95.01327 92.95829 90.96206 91.96742 90.12840 93.46854 93.47244
## 2004 95.49327 98.81829 97.24206 97.90742 97.41840 96.70854 97.82244
## 2005 97.29327 96.53829 97.58206 98.49742 96.04840 101.19854 100.87244
## 2006 104.65327 105.81829 108.07206 109.85742 107.80840 111.69854 109.52244
## 2007 109.48327 110.05829 117.65206 111.03742 113.17840 115.33854 110.42244
## 2008 114.91327 111.72829 113.32206 114.96742 107.94840 110.08854 105.16244
## 2009 82.89327 81.09829 78.33206 80.43742 78.90840 77.82854 81.38244
## 2010 84.72327 84.48829 86.34206 90.79742 89.73840 95.42854 90.96244
## 2011 97.03327 93.74829 95.35206 94.90742 97.13840 94.45854 94.34244
## Aug Sep Oct Nov Dec
## 1996 79.34631 78.14206 80.25897 80.52212 78.34012
## 1997 86.74631 83.82206 86.58897 87.35212 81.52012
## 1998 83.32631 86.09206 83.09897 83.14212 81.97012
## 1999 96.09631 96.32206 98.06897 95.75212 100.49012
## 2000 113.10631 116.29206 113.30897 113.23212 117.30012
## 2001 96.68631 94.82206 93.34897 97.61212 98.52012
## 2002 92.60631 92.24206 93.63897 92.45212 89.67012
## 2003 93.88631 97.14206 96.89897 97.33212 97.71012
## 2004 96.48631 98.19206 96.21897 96.04212 102.27012
## 2005 100.17631 103.49206 102.05897 103.33212 102.93012
## 2006 107.37631 114.49206 114.18897 109.89212 108.90012
## 2007 110.85631 114.58206 113.85897 112.01212 115.83012
## 2008 106.09631 102.05206 100.99897 96.30212 89.70012
## 2009 82.11631 80.56206 81.51897 84.16212 90.97012
## 2010 94.53631 88.32206 91.09897 97.22212 89.81012
## 2011 93.31631 86.64206 88.73897 88.21212plot(eeadj)
The time plot shows some sudden changes, particularly the big drop in 2008/2009. These changes are due to the global economic environment. Otherwise there is nothing unusual about the time plot and there appears to be no need to do any data adjustments.
There is no evidence of changing variance, so we will not do a Box-Cox transformation.
The data are clearly non-stationary as the series wanders up and down for long periods. Consequently, we will take a first difference of the data. The differenced data are shown in Figure 8.12. These look stationary, and so we will not consider further differences.
tsdisplay(diff(eeadj),main="")
The PACF shown in Figure 8.12 is suggestive of an AR(3) model. So an initial candidate model is an ARIMA(3,1,0). There are no other obvious candidate models.
We fit an ARIMA(3,1,0) model along with variations including ARIMA(4,1,0), ARIMA(2,1,0), ARIMA(3,1,1), etc. Of these, the ARIMA(3,1,1) has a slightly smaller AICc value.
fit <- Arima(eeadj, order=c(3,1,1))
summary(fit)## Series: eeadj
## ARIMA(3,1,1)
##
## Coefficients:
## ar1 ar2 ar3 ma1
## 0.0519 0.1191 0.3730 -0.4542
## s.e. 0.1840 0.0888 0.0679 0.1993
##
## sigma^2 estimated as 9.737: log likelihood=-484.08
## AIC=978.17 AICc=978.49 BIC=994.4
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.001227744 3.079373 2.389267 -0.04290849 2.517748 0.2913919
## ACF1
## Training set 0.008928479- The ACF plot of the residuals from the ARIMA(3,1,1) model shows all correlations within the threshold limits indicating that the residuals are behaving like white noise. A portmanteau test returns a large p-value, also suggesting the residuals are white noise.
Acf(residuals(fit))
Box.test(residuals(fit), lag=24, fitdf=4, type="Ljung")##
## Box-Ljung test
##
## data: residuals(fit)
## X-squared = 20.496, df = 20, p-value = 0.4273High p-value suggests residuals are white noise
Now I will try the auto.arima on the original data
model2=auto.arima(elecequip)
summary(model2)## Series: elecequip
## ARIMA(4,0,1)(0,1,1)[12]
##
## Coefficients:
## ar1 ar2 ar3 ar4 ma1 sma1
## 1.1661 -0.0149 0.2053 -0.3825 -0.6083 -0.7975
## s.e. 0.1322 0.1348 0.1210 0.0718 0.1376 0.0738
##
## sigma^2 estimated as 11.09: log likelihood=-473.73
## AIC=961.46 AICc=962.12 BIC=983.77
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set 0.170688 3.169897 2.373325 0.09769198 2.471343 0.2894476
## ACF1
## Training set 0.008846644# this is not important, but be familiar with it just in case
plot(model2)
plot(model2$residuals)
Acf(residuals(model2))
Box.test(residuals(model2), lag=24, fitdf=4, type="Ljung")##
## Box-Ljung test
##
## data: residuals(model2)
## X-squared = 16.881, df = 20, p-value = 0.6607Now let’s look at forecasts:
plot(forecast(fit))
plot(forecast(model2))
5 Week5 Overview
The “frequency” is the number of observations per “cycle” (normally a year, but sometimes a week, a day, an hour, etc). This is the opposite of the definition of frequency in physics, or in Fourier analysis, where “period” is the length of the cycle, and “frequency” is the inverse of period. When using the ts() function in R, the following choices should be used.
Data frequency
Annual 1
Quarterly 4
Monthly 12
Weekly 52
Actually, there are not 52 weeks in a year, but 365.25/7 = 52.18 on average. But most functions which use ts objects require integer frequency.
Once the frequency of observations is smaller than a week, then there is usually more than one way of handling the frequency. For example, data observed every minute might have an hourly seasonality (frequency=60), a daily seasonality (frequency=24x60=1440), a weekly seasonality (frequency=24x60x7=10080) and an annual seasonality (frequency=24x60x365.25=525960). If you want to use a ts object, then you need to decide which of these is the most important.
An alternative is to use a msts object (defined in the forecast package) which handles multiple seasonality time series. Then you can specify all the frequencies that might be relevant. It is also flexible enough to handle non-integer frequencies.
You won’t necessarily want to include all of these frequencies — just the ones that are likely to be present in the data. As you have only 180 days of data, you can probably ignore the annual seasonality. If the data are measurements of a natural phenomenon (e.g., temperature), you might also be able to ignore the weekly seasonality.
With multiple seasonalities, you could use a TBATS model, or Fourier terms in a regression or ARIMA model. The fourier function from the forecast package will handle msts objects.
5.1 weekly forecasting
This is another situation where Fourier terms are useful for handling the seasonality. Not only is the seasonal period rather long, it is non-integer (averaging 365.25/7 = 52.18). So ARIMA and ETS models do not tend to give good results, even with a period of 52 as an approximation.
library(forecast)
gas <- ts(read.csv("https://robjhyndman.com/data/gasoline.csv", header=FALSE)[,1],
freq=365.25/7, start=1991+31/365.25)
bestfit <- list(aicc=Inf)
for(i in 1:25)
{
fit <- auto.arima(gas, xreg=fourier(gas, K=i), seasonal=FALSE)
if(fit$aicc < bestfit$aicc)
bestfit <- fit
else break;
}
fc <- forecast(bestfit, xreg=fourier(gas, K=12, h=104))
plot(fc)
or An alternative approach is the TBATS model introduced by De Livera et al (JASA, 2011). This uses a state space model that is a generalization of those underpinning exponential smoothing. It also allows for automatic Box-Cox transformation and ARMA errors. The modelling algorithm is entirely automated:
gastbats <- tbats(gas)
fc2 <- forecast(gastbats, h=104)
plot(fc2, ylab="thousands of barrels per day")
Here the fitted model is given at the top of the plot as TBATS(0.999, {2,2}, 1, {<52.18,8>}). That is, a Box-Cox transformation of 0.999 (essentially doing nothing), ARMA(2,2) errors, a damping parameter of 1 (doing nothing) and 8 Fourier pairs with period m = 52.18 .
5.2 Time Series validation:
5.2.1 k-fold Cross-Validation applied to an autoregressive model
library(forecast)
library(fpp)
modelcv <- CVar(lynx, k=5, lambda=0.15)
print(modelcv)## Series: lynx
## Call: CVar(y = lynx, k = 5, lambda = 0.15)
##
## 5-fold cross-validation
## Mean SD
## ME -9.2398304 112.7521276
## RMSE 853.3043240 219.8854819
## MAE 564.9275194 153.7261921
## MPE -7.3619245 4.2851572
## MAPE 46.1519401 6.9976745
## ACF1 0.0839324 0.1992721
## Theil's U 0.7720786 0.1151928Time Series Cross-Validation:
differnt ways to think about this:
Fixing the start,
Rolling Windows
Future horizon
Quatrterly moving, etc.
library(fpp)
e <- tsCV(dj, rwf, drift=TRUE, h=1)
sqrt(mean(e^2, na.rm=TRUE))## [1] 22.68249## [1] 22.68249
sqrt(mean(residuals(rwf(dj, drift=TRUE))^2, na.rm=TRUE))## [1] 22.49681## [1] 22.49681library(fpp) # To load the data set a10
plot(a10, ylab="$ million", xlab="Year", main="Antidiabetic drug sales")
plot(log(a10), ylab="", xlab="Year", main="Log Antidiabetic drug sales")
k <- 60 # minimum data length for fitting a model
n <- length(a10)
mae1 <- mae2 <- mae3 <- mae4<- matrix(NA,12,12)
st <- tsp(a10)[1]+(k-1)/12
for(i in 1:12)
{
xshort <- window(a10, end=st + (i-1))
xnext <- window(a10, start=st + (i-1) + 1/12, end=st + i)
fit1 <- tslm(xshort ~ trend + season, lambda=0)
fcast1 <- forecast(fit1, h=12)
fit2 <- Arima(xshort, order=c(3,0,1), seasonal=list(order=c(0,1,1), period=12),
include.drift=TRUE, lambda=0, method="ML")
fcast2 <- forecast(fit2, h=12)
fit3 <- ets(xshort,model="MMM",damped=TRUE)
fcast3 <- forecast(fit3, h=12)
fit4 <- tbats(xshort)
fcast4 <- forecast(fit4, h=12)
mae1[i,] <- abs(fcast1[['mean']]-xnext)
mae2[i,] <- abs(fcast2[['mean']]-xnext)
mae3[i,] <- abs(fcast3[['mean']]-xnext)
mae4[i,] <- abs(fcast4[['mean']]-xnext)
}
plot(1:12, colMeans(mae1), type="l", col=2, xlab="horizon", ylab="MAE",
ylim=c(0.35,1.5))
lines(1:12, colMeans(mae2), type="l",col=3)
lines(1:12, colMeans(mae3), type="l",col=4)
lines(1:12, colMeans(mae4), type="l",col=5)
legend("topleft",legend=c("LM","ARIMA","ETS","Tbats"),col=2:5,lty=1)
6 Week6 Overview
6.1 Agenda
Regression with ARMA errors
ARMAX models
Long memory ARIMA models
Forecasting non-zero mean processes
6.2 Dynamic Regression and Regression with ARMA errors
In regression error should not be autocorrelated (it has to be random, right?)
What If \(e_t\) is not random and has autocorrlation. In that case we can re write the regression with ARIMA error modeling to capture this behaviour. By doing so the new error term will become White noise.
When we estimate the parameters from the model, we need to minimize the sum of squared \(e_t\) values. If, instead, we minimized the sum of squared \(n_t\) values (which is what would happen if we estimated the regression model ignoring the autocorrelations in the errors), then several problems arise.
The estimated coefficients β0,…,βk are no longer the best estimates as some information has been ignored in the calculation;
Statistical tests associated with the model (e.g., t-tests on the coefficients) are incorrect.
The AIC of the fitted models are not a good guide as to which is the best model for forecasting.
6.2.1 Spurious Regression
In most cases, the pp-values associated with the coefficients will be too small, and so some predictor variables appear to be important when they are not. This is known as “spurious regression”.
Note:An important consideration in estimating a regression with ARMA errors is that all variables in the model must first be stationary. So we first have to check that ytyt and all the predictors (x1,t,…,xk,t) appear to be stationary. If we estimate the model while any of these are non-stationary, the estimated coefficients can be incorrect.One exception to this is the case where non-stationary variables are co-integrated. If there exists a linear combination between the non-stationary ytyt and predictors that is stationary, then the estimated coefficients are correct.
6.2.2 Regression with ARIMA errors modeling :
Check that the forecast variable and all predictors are stationary. If not, apply differencing until all variables are stationary. Where appropriate, use the same differencing for all variables to preserve interpretability.
Fit the regression model with AR(2) errors for non-seasonal data or ARIMA(2,0,0)(1,0,0)mm errors for seasonal data.
Calculate the errors (ntnt) from the fitted regression model and identify an appropriate ARMA model for them.
Re-fit the entire model using the new ARMA model for the errors. Check that the etet series looks like white noise.
library(fpp)
plot(usconsumption, xlab="Year",
main="Quarterly changes in US consumption and personal income")
head(usconsumption)## consumption income
## 1970 Q1 0.6122769 0.496540
## 1970 Q2 0.4549298 1.736460
## 1970 Q3 0.8746730 1.344881
## 1970 Q4 -0.2725144 -0.328146
## 1971 Q1 1.8921870 1.965432
## 1971 Q2 0.9133782 1.490757We would like to forecast changes in consumption is based on changes in income.
cor(usconsumption[,1], usconsumption[,2])## [1] 0.4320562fit <- Arima(usconsumption[,1], xreg=usconsumption[,2],
order=c(2,0,0))
tsdisplay(arima.errors(fit), main="ARIMA errors")## Deprecated, use residuals.Arima(object, type='regression') instead
(fit2 <- Arima(usconsumption[,1], xreg=usconsumption[,2],
order=c(1,0,2)))## Series: usconsumption[, 1]
## Regression with ARIMA(1,0,2) errors
##
## Coefficients:
## ar1 ma1 ma2 intercept usconsumption[, 2]
## 0.6516 -0.5440 0.2187 0.5750 0.2420
## s.e. 0.1468 0.1576 0.0790 0.0951 0.0513
##
## sigma^2 estimated as 0.3502: log likelihood=-144.27
## AIC=300.54 AICc=301.08 BIC=319.14Now you can show residuals are White Noise
Box.test(residuals(fit2),fitdf=5,lag=10,type="Ljung")##
## Box-Ljung test
##
## data: residuals(fit2)
## X-squared = 4.5948, df = 5, p-value = 0.4673Forecasts are, of course, only possible if we have future values of changes in personal disposable income. We assume for the next 8 quarters, the percentage change in personal disposable income is equal to the mean percentage change from the last forty years.
fcast <- forecast(fit2,xreg=rep(mean(usconsumption[,2]),8), h=8)
plot(fcast,
main="Forecasts from regression with ARIMA(1,0,2) errors")
we can also do this with auto.arima()
fit <- auto.arima(usconsumption[,1], xreg=usconsumption[,2])To forecast a regression model with ARIMA errors, we need to forecast the regression part of the model and the ARIMA part of the model, and combine the results.
example2:
plot(insurance, main="Insurance advertising and quotations", xlab="Year")
# Lagged predictors. Test 0, 1, 2 or 3 lags.
Advert <- cbind(insurance[,2],
c(NA,insurance[1:39,2]),
c(NA,NA,insurance[1:38,2]),
c(NA,NA,NA,insurance[1:37,2]))
colnames(Advert) <- paste("AdLag",0:3,sep="")
# Choose optimal lag length for advertising based on AIC
# Restrict data so models use same fitting period
fit1 <- auto.arima(insurance[4:40,1], xreg=Advert[4:40,1], d=0)
fit2 <- auto.arima(insurance[4:40,1], xreg=Advert[4:40,1:2], d=0)
fit3 <- auto.arima(insurance[4:40,1], xreg=Advert[4:40,1:3], d=0)
fit4 <- auto.arima(insurance[4:40,1], xreg=Advert[4:40,1:4], d=0)
# Best model fitted to all data (based on AICc)
# Refit using all data
fit <- auto.arima(insurance[,1], xreg=Advert[,1:2], d=0)
fit## Series: insurance[, 1]
## Regression with ARIMA(3,0,0) errors
##
## Coefficients:
## ar1 ar2 ar3 intercept AdLag0 AdLag1
## 1.4117 -0.9317 0.3591 2.0393 1.2564 0.1625
## s.e. 0.1698 0.2545 0.1592 0.9931 0.0667 0.0591
##
## sigma^2 estimated as 0.2165: log likelihood=-23.89
## AIC=61.78 AICc=65.28 BIC=73.66.3 ARIMAX models vs. Regression with ARMA errors
They are differnt, both of these models can be considered as special cases of transfer function models, popularized by Box and Jenkins.
An ARMAX model simply adds in the covariate on the right hand side, which make the coeffiecnets not really unterpretable
6.4 Autoregressive Fractionally Integrated Moving Average (ARFIMA):
The long lasting autocorrelation function of the data shows the presence of long memory structure, and the Hurst exponent test could be used to test for presence of long memory structure. The Geweke and Porter-Hudak (GPH) method of estimation could be used to obtain the long memory parameter d of the ARFIMA model. In ARFIMA d is less than. The ARFIMA model searches for a non-integer parameter, d, to differentiate the data to capture long memory.
This model allows for the series to be fractionally integrated, generalizing the ARIMA model’s integer order of integration to allow the d parameter to take on fractional values, \(−0.5 < d < 0.5\).
In estimating an ARIMA model, the researcher chooses the integer order of differencing d to ensure that the resulting series stationary process. As unit root tests often lack the power to distinguish between a truly nonstationary (I(1)) series and a stationary series embodying a structural break or shift, time series are often first-differenced if they do not receive a clean bill of health from unit root testing. Many time series exhibit too much long-range dependence to be classified as I(0) but are not I(1). The ARFIMA model is designed to represent these series.
7 Week7 Overview
7.1 Agenda
7.2 stat Space models (dynamic linear models, DLM)
State-space models represent the role of historyin a finite-dimensional vector (in other words they don’t look at the whole process), the state. Very Markovian. Autoregressions are a special case in which the state vector — the previous p observations — is observed.
Markov Chain reminder:
\(p(a,b,c,d)=p(a)*p(b|a)*p(c|b)*p(d|c)\)
A linear-Gaussian state space model for an m−dimensional time series yt consists of a measurement equation relating the observed data to an p−dimensional state vector θt, and a Markovian transition equation that describes the evolution of the state vector over time
7.3 Vector autoregressions (VAR) and VARIMA
This is important where all variables affect each other. In VAR all variables are treated symmetrically. They are all modelled as if they influence each other equally. In other words all variables are now treated as “endogenous”(Endogenous variables have values that are determined by other variables in the system)
library(vars)## Warning: package 'vars' was built under R version 3.4.4## Loading required package: MASS##
## Attaching package: 'MASS'## The following objects are masked from 'package:fma':
##
## cement, housing, petrol## Loading required package: strucchange## Warning: package 'strucchange' was built under R version 3.4.4## Loading required package: sandwich## Loading required package: urca## Warning: package 'urca' was built under R version 3.4.4library(fpp)
VARselect(usconsumption, lag.max=8, type="const")$selection## AIC(n) HQ(n) SC(n) FPE(n)
## 5 1 1 5var <- VAR(usconsumption, p=3, type="const")
serial.test(var, lags.pt=10, type="PT.asymptotic")##
## Portmanteau Test (asymptotic)
##
## data: Residuals of VAR object var
## Chi-squared = 33.384, df = 28, p-value = 0.2219summary(var)##
## VAR Estimation Results:
## =========================
## Endogenous variables: consumption, income
## Deterministic variables: const
## Sample size: 161
## Log Likelihood: -338.797
## Roots of the characteristic polynomial:
## 0.7666 0.5529 0.524 0.524 0.3181 0.3181
## Call:
## VAR(y = usconsumption, p = 3, type = "const")
##
##
## Estimation results for equation consumption:
## ============================================
## consumption = consumption.l1 + income.l1 + consumption.l2 + income.l2 + consumption.l3 + income.l3 + const
##
## Estimate Std. Error t value Pr(>|t|)
## consumption.l1 0.22280 0.08580 2.597 0.010326 *
## income.l1 0.04037 0.06230 0.648 0.518003
## consumption.l2 0.20142 0.09000 2.238 0.026650 *
## income.l2 -0.09830 0.06411 -1.533 0.127267
## consumption.l3 0.23512 0.08824 2.665 0.008530 **
## income.l3 -0.02416 0.06139 -0.394 0.694427
## const 0.31972 0.09119 3.506 0.000596 ***
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Residual standard error: 0.6304 on 154 degrees of freedom
## Multiple R-Squared: 0.2183, Adjusted R-squared: 0.1878
## F-statistic: 7.166 on 6 and 154 DF, p-value: 9.384e-07
##
##
## Estimation results for equation income:
## =======================================
## income = consumption.l1 + income.l1 + consumption.l2 + income.l2 + consumption.l3 + income.l3 + const
##
## Estimate Std. Error t value Pr(>|t|)
## consumption.l1 0.48705 0.11637 4.186 4.77e-05 ***
## income.l1 -0.24881 0.08450 -2.945 0.003736 **
## consumption.l2 0.03222 0.12206 0.264 0.792135
## income.l2 -0.11112 0.08695 -1.278 0.203170
## consumption.l3 0.40297 0.11967 3.367 0.000959 ***
## income.l3 -0.09150 0.08326 -1.099 0.273484
## const 0.36280 0.12368 2.933 0.003865 **
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Residual standard error: 0.855 on 154 degrees of freedom
## Multiple R-Squared: 0.2112, Adjusted R-squared: 0.1805
## F-statistic: 6.873 on 6 and 154 DF, p-value: 1.758e-06
##
##
##
## Covariance matrix of residuals:
## consumption income
## consumption 0.3974 0.1961
## income 0.1961 0.7310
##
## Correlation matrix of residuals:
## consumption income
## consumption 1.0000 0.3639
## income 0.3639 1.0000fcst <- forecast(var)
plot(fcst, xlab="Year")
8 FAQ
Why \(r^2\) is not valid for non-linear regression?
http://statisticsbyjim.com/regression/r-squared-invalid-nonlinear-regression/
What is the mean and variance of White Noise?
White noise has zero mean and finite variance.
How to deal with Multicollinearity
What is innovation? The innovation is the difference between the observed value of a variable at time t and the optimal forecast of that value based on information available prior to time t. If the forecasting method is working correctly, successive innovations are uncorrelated with each other (i.e. white noise).
Thus it can be said that the innovation time series is obtained from the measurement time series by a process of ‘whitening’, or removing the predictable component.
what is drift?
In probability theory, stochastic drift is the change of the average value of a stochastic (random) process.
Reference: