Question 1:
library("Quandl")
## Loading required package: xts
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
library("tseries")
library("forecast")
## Loading required package: timeDate
## This is forecast 6.2
cpil <-Quandl("FRED/CPILFENS", type ="ts")
str(cpil)
## Time-Series [1:707] from 1957 to 2016: 28.5 28.5 28.7 28.8 28.8 28.9 28.9 29 29.1 29.2 ...
plot(cpil)
lncpil<-log(cpil)
inflation <-diff(lncpil,12)
plot(inflation)
inflationall <-inflation
inf1 <-window(inflationall, end=c(2014,12))
inf2 <-window(inflationall, start=c(2015,1))
library(tseries)
adf.test(inf1)
##
## Augmented Dickey-Fuller Test
##
## data: inf1
## Dickey-Fuller = -3.1898, Lag order = 8, p-value = 0.08969
## alternative hypothesis: stationary
dinf1<-diff(inf1)
plot(dinf1)
adf.test(dinf1)
## Warning in adf.test(dinf1): p-value smaller than printed p-value
##
## Augmented Dickey-Fuller Test
##
## data: dinf1
## Dickey-Fuller = -6.211, Lag order = 8, p-value = 0.01
## alternative hypothesis: stationary
library(tseries)
kpss.test(inf1,null ="Trend")
## Warning in kpss.test(inf1, null = "Trend"): p-value smaller than printed p-
## value
##
## KPSS Test for Trend Stationarity
##
## data: inf1
## KPSS Trend = 1.3584, Truncation lag parameter = 6, p-value = 0.01
kpss.test(dinf1, null ="Trend")
## Warning in kpss.test(dinf1, null = "Trend"): p-value greater than printed
## p-value
##
## KPSS Test for Trend Stationarity
##
## data: dinf1
## KPSS Trend = 0.046471, Truncation lag parameter = 6, p-value = 0.1
dinf1 <-diff(inf1)
plot(dinf1)
adf.test(dinf1)
## Warning in adf.test(dinf1): p-value smaller than printed p-value
##
## Augmented Dickey-Fuller Test
##
## data: dinf1
## Dickey-Fuller = -6.211, Lag order = 8, p-value = 0.01
## alternative hypothesis: stationary
par(mfrow = c(2,2))
acf(as.data.frame(dinf1),type= 'correlation',lag = 36,main="Sample ACF")
acf(as.data.frame(dinf1),type = 'partial',lag =36, main=" Simple PACF")
dinf <-diff(inflation)
am <-arima(dinf,order=c(3,1,3), seasonal=list(order=c(0,0,2),period=12))
am
##
## Call:
## arima(x = dinf, order = c(3, 1, 3), seasonal = list(order = c(0, 0, 2), period = 12))
##
## Coefficients:
## ar1 ar2 ar3 ma1 ma2 ma3 sma1 sma2
## -0.7591 -0.7009 0.2062 0.0936 0.1352 -0.7966 -0.8058 -0.0283
## s.e. 0.0606 0.0744 0.0508 0.0461 0.0521 0.0414 0.0408 0.0411
##
## sigma^2 estimated as 3.116e-06: log likelihood = 3401.97, aic = -6785.95
tsdiag(am,gof.lag=24)
am.LB <-Box.test(am$residuals, lag=24, type="Ljung")
am.LB
##
## Box-Ljung test
##
## data: am$residuals
## X-squared = 28.839, df = 24, p-value = 0.2262
M4 <-arima(dinf,order=c(1,0,4),seasonal=list(order=c(2,0,2),period=12))
M4
##
## Call:
## arima(x = dinf, order = c(1, 0, 4), seasonal = list(order = c(2, 0, 2), period = 12))
##
## Coefficients:
## ar1 ma1 ma2 ma3 ma4 sar1 sar2 sma1
## 0.9674 -0.6606 -0.0098 -0.1021 -0.0078 -0.3516 0.0226 -0.4440
## s.e. 0.0141 0.0408 0.0489 0.0479 0.0404 3.1234 0.1408 3.1225
## sma2 intercept
## -0.3274 0e+00
## s.e. 2.6074 1e-04
##
## sigma^2 estimated as 3.109e-06: log likelihood = 3409.22, aic = -6796.43
tsdiag(M4,gof.lag=24)
M4.LB <-Box.test(M4$residuals, lag=24, type="Ljung")
par(mfrow=c(1,2), cex=0.75)
am.forecast <- forecast(am, h=24)
plot(am.forecast, xlim=c(2011,2016))