Homework 5, Problem 1
Andrii Baryshpolets
30 April 2017
Intervention analysis: Retail and Food Services Sales
options(allow_html_dependencies=TRUE)
library(Quandl)## Warning: package 'Quandl' was built under R version 3.3.3## Loading required package: xts## Loading required package: zoo##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numericlibrary(zoo)
library(xts)
library(dygraphs)
library(knitr)
library(forecast)## Loading required package: timeDate## This is forecast 7.3library(urca)
Quandl.api_key('rxVQhZ8_nxxeo2yy4Uz4')
yall <- Quandl("FRED/RSAFSNA", type="ts") (1)
lyall <- log(yall)
par(mfrow=c(2,1), cex=0.85)
plot(yall, xlab=" Years 1992-2017", ylab="", main="Retail and Food Services Sales" )
plot(lyall, xlab="Years 1992-2017", ylab="", main="Log of Retail and Food Services Sales" )
y <- window(yall, end=c(2008 + 5/12))
ly <- log(y)
dly1 <- diff(ly)
dly12 <- diff(ly, lag=12)
d2ly1_12 <- diff(diff(ly), lag=12)
dygraph(dly1)maxlag <- 60
par(mfrow=c(3,4))
plot(ly, ylab='', xlab='', main='log(RFSS)')
plot(dly1, ylab='', xlab='year', main=expression(paste(Delta, "log(RFSS)")))
plot(dly12, ylab='', xlab='year', main=expression(paste(Delta[12], "log(RFSS)")))
plot(d2ly1_12, ylab='', xlab='year', main=expression(paste(Delta,Delta[12], "log(RFSS)")))
Acf(ly, type="correlation", lag.max=maxlag, main=expression(paste("ACF for log(RFSS)")))
Acf(dly1, type="correlation", lag.max=maxlag, main=expression(paste("ACF for ",Delta, "log(RFSS)")))
Acf(dly12, type="correlation", lag.max=maxlag, main=expression(paste("ACF for ",Delta[12], "log(RFSS)")))
Acf(d2ly1_12, type="correlation", lag.max=maxlag, main=expression(paste("ACF for ",Delta,Delta[12], "log(RFSS)")))
Acf(ly, type="partial", lag.max=maxlag, main=expression(paste("PACF for log(RFSS)")))
Acf(dly1, type="partial", lag.max=maxlag, main=expression(paste("PACF for ",Delta, "log(RFSS)")))
Acf(dly12, type="partial", lag.max=maxlag, main=expression(paste("PACF for ",Delta[12], "log(RFSS)")))
Acf(d2ly1_12, type="partial", lag.max=maxlag, main=expression(paste("PACF for ",Delta,Delta[12], "log(RFSS)")))
We use autoarima function to find the best specification for our model.
M1.bic <- auto.arima(y, ic="bic", seasonal=TRUE, stationary=FALSE, stepwise=FALSE)
M1.bic## Series: y
## ARIMA(3,0,0)(0,1,1)[12] with drift
##
## Coefficients:
## ar1 ar2 ar3 sma1 drift
## 0.1438 0.1958 0.4945 -0.5061 1081.1111
## s.e. 0.0641 0.0639 0.0652 0.0683 91.2315
##
## sigma^2 estimated as 25122355: log likelihood=-1848.25
## AIC=3708.5 AICc=3708.97 BIC=3727.85M1.aic <- auto.arima(y, ic="aic", seasonal=TRUE, stationary=FALSE, stepwise=FALSE)
M1.aic## Series: y
## ARIMA(3,0,0)(0,1,2)[12] with drift
##
## Coefficients:
## ar1 ar2 ar3 sma1 sma2 drift
## 0.1812 0.1816 0.4319 -0.3733 -0.1691 1087.890
## s.e. 0.0681 0.0672 0.0735 0.0857 0.0828 71.259
##
## sigma^2 estimated as 24778313: log likelihood=-1846.33
## AIC=3706.66 AICc=3707.29 BIC=3729.24M1.aicc <- auto.arima(y, ic="aicc", seasonal=TRUE, stationary=FALSE, stepwise=FALSE)
M1.aicc## Series: y
## ARIMA(3,0,0)(0,1,2)[12] with drift
##
## Coefficients:
## ar1 ar2 ar3 sma1 sma2 drift
## 0.1812 0.1816 0.4319 -0.3733 -0.1691 1087.890
## s.e. 0.0681 0.0672 0.0735 0.0857 0.0828 71.259
##
## sigma^2 estimated as 24778313: log likelihood=-1846.33
## AIC=3706.66 AICc=3707.29 BIC=3729.24Two of three criterias suggest to use ARIMA(3,0,0)(0,1,2)[12].
m2 <- Arima(y, order=c(3,0,0), seasonal=list(order=c(0,1,2), period=12))
m2## Series: y
## ARIMA(3,0,0)(0,1,2)[12]
##
## Coefficients:
## ar1 ar2 ar3 sma1 sma2
## 0.2361 0.2480 0.5070 -0.4047 -0.1285
## s.e. 0.0676 0.0634 0.0681 0.0876 0.0845
##
## sigma^2 estimated as 25729530: log likelihood=-1851.26
## AIC=3714.51 AICc=3714.98 BIC=3733.86Forecast Using Pre-Intervention Model
m2.f <- forecast(m2, h=78)
y.1214 <- window(yall, end=c(2014,12))
m2.f.err <- y.1214 - m2.f$meanWe plot forecast, forecast error ACF and PACF for forecast errors
plot(m2.f, xlim=c(2007, 2015), ylim=c(300000,560000))
lines(index(y.1214), y.1214)
plot(m2.f.err, ylim=c(-55000,20000), main="Multistep Forecast Error")
abline(h=0, lty="dashed")
par(mfrow=c(2,1), mar=c(2,4,2,2), cex=0.9)
Acf(m2.f.err, type="correlation", lag.max=48, ylab="ACF", main="")
Acf(m2.f.err, type="partial", lag.max=48, ylab="PACF", main="")
(2)
library(TSA)## Warning: package 'TSA' was built under R version 3.3.3## Loading required package: leaps## Warning: package 'leaps' was built under R version 3.3.3## Loading required package: locfit## Warning: package 'locfit' was built under R version 3.3.3## locfit 1.5-9.1 2013-03-22## Loading required package: mgcv## Loading required package: nlme##
## Attaching package: 'nlme'## The following object is masked from 'package:forecast':
##
## getResponse## This is mgcv 1.8-15. For overview type 'help("mgcv-package")'.## Loading required package: tseries##
## Attaching package: 'TSA'## The following objects are masked from 'package:forecast':
##
## fitted.Arima, plot.Arima## The following objects are masked from 'package:timeDate':
##
## kurtosis, skewness## The following objects are masked from 'package:stats':
##
## acf, arima## The following object is masked from 'package:utils':
##
## tarlibrary(cpm)
y.1214 <- window(yall, end=c(2014,12))
Pjun8 <- 1*(index(y.1214)==2008 + 5/12)
m3 <- arimax(log(y.1214), order=c(3,0,0), seasonal=list(order=c(0,1,2), period=12), xtransf=data.frame(Pjun8), transfer=list(c(1,0)), method="ML")
m3##
## Call:
## arimax(x = log(y.1214), order = c(3, 0, 0), seasonal = list(order = c(0, 1,
## 2), period = 12), method = "ML", xtransf = data.frame(Pjun8), transfer = list(c(1,
## 0)))
##
## Coefficients:
## ar1 ar2 ar3 sma1 sma2 Pjun8-AR1 Pjun8-MA0
## 0.4150 0.2358 0.3466 -0.5618 -0.1959 0.6953 0.0232
## s.e. 0.0613 0.0618 0.0623 0.0660 0.0649 0.1368 0.0154
##
## sigma^2 estimated as 0.0003597: log likelihood = 666.52, aic = -1319.04m4 <- arimax(log(y.1214), order=c(3,0,0), seasonal=list(order=c(0,1,2), period=12),
xtransf=data.frame(Pjun8,Pjun8), transfer=list(c(0,0),c(1,0)), method="ML")
m4##
## Call:
## arimax(x = log(y.1214), order = c(3, 0, 0), seasonal = list(order = c(0, 1,
## 2), period = 12), method = "ML", xtransf = data.frame(Pjun8, Pjun8), transfer = list(c(0,
## 0), c(1, 0)))
##
## Coefficients:
## ar1 ar2 ar3 sma1 sma2 Pjun8-MA0 Pjun8.1-AR1
## 0.4389 0.2152 0.3433 -0.5482 -0.2123 0.8258 -0.0332
## s.e. 0.0602 0.0633 0.0619 0.0652 0.0637 2.3233 0.0909
## Pjun8.1-MA0
## -0.8241
## s.e. 2.3224
##
## sigma^2 estimated as 0.000358: log likelihood = 667.12, aic = -1318.25tsdiag(m3, gof.lag=24)
tsdiag(m4, gof.lag=24)
(3)
library(TSA)
library(cpm)
hmax <- 24
Pjun8 <- 1*(index(y.1214)==2008 + 5/12)
Pjun8x <- c(Pjun8, rep(0,hmax))
tf4 <- cbind(m4$coef["Pjun8-MA0"]*Pjun8x, m4$coef["Pjun8.1-MA0"]*filter(Pjun8x, filter=m4$coef["Pjun8.1-AR1"], method='recursive', sides=1))
m4x <- Arima(log(y.1214), order=c(3,0,0), seasonal=list(order=c(0,1,2), period=12), xreg=tf4[1:length(y.1214),])
m4x## Series: log(y.1214)
## ARIMA(3,0,0)(0,1,2)[12]
##
## Coefficients:
## ar1 ar2 ar3 sma1 sma2
## 0.4383 0.2148 0.3443 -0.5495 -0.2113
## s.e. 0.0602 0.0632 0.0618 0.0652 0.0638
## m4$coef["Pjun8-MA0"] * Pjun8x
## 1.0220
## s.e. 0.5556
## m4$coef["Pjun8.1-MA0"] * filter(Pjun8x, filter = m4$coef["Pjun8.1-AR1"],
## 1.0219
## s.e. 0.5526
##
## sigma^2 estimated as 0.0003745: log likelihood=667.12
## AIC=-1318.25 AICc=-1317.69 BIC=-1289.64m4x.f.h <- forecast(m4x, h=hmax, xreg=tf4[(length(y.1214)+1):(length(y.1214)+hmax),])Plot a multistep forecast for the model with transfer function.
plot(m4x.f.h, xlim=c(2014, 2017), ylim=c(12.8,13.3))
lines(log(yall), col="black", lwd=2, xlim=c(2012, 2017), ylim=c(12.5,14))