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)
#(1) data identification
unr <- Quandl("FRED/UNRATE", type="zoo")
summary(unr)
## Index unr
## Min. :1948 Min. : 2.500
## 1st Qu.:1965 1st Qu.: 4.700
## Median :1982 Median : 5.600
## Mean :1982 Mean : 5.824
## 3rd Qu.:1999 3rd Qu.: 6.900
## Max. :2016 Max. :10.800
str(unr)
## 'zooreg' series from Jan 1948 to Jan 2016
## Data: num [1:817] 3.4 3.8 4 3.9 3.5 3.6 3.6 3.9 3.8 3.7 ...
## Index: Class 'yearmon' num [1:817] 1948 1948 1948 1948 1948 ...
## Frequency: 12
head(unr)
## Jan 1948 Feb 1948 Mar 1948 Apr 1948 May 1948 Jun 1948
## 3.4 3.8 4.0 3.9 3.5 3.6
tail(unr)
## Aug 2015 Sep 2015 Oct 2015 Nov 2015 Dec 2015 Jan 2016
## 5.1 5.1 5.0 5.0 5.0 4.9
y <- unr
d.y <-diff(y)
plot(y)
plot(d.y)
adf.test(y)
## Warning in adf.test(y): p-value smaller than printed p-value
##
## Augmented Dickey-Fuller Test
##
## data: y
## Dickey-Fuller = -4.0267, Lag order = 9, p-value = 0.01
## alternative hypothesis: stationary
kpss.test(y)
## Warning in kpss.test(y): p-value smaller than printed p-value
##
## KPSS Test for Level Stationarity
##
## data: y
## KPSS Level = 1.9577, Truncation lag parameter = 6, p-value = 0.01
adf.test(d.y)
## Warning in adf.test(d.y): p-value smaller than printed p-value
##
## Augmented Dickey-Fuller Test
##
## data: d.y
## Dickey-Fuller = -8.2434, Lag order = 9, p-value = 0.01
## alternative hypothesis: stationary
kpss.test(d.y)
## Warning in kpss.test(d.y): p-value greater than printed p-value
##
## KPSS Test for Level Stationarity
##
## data: d.y
## KPSS Level = 0.064885, Truncation lag parameter = 6, p-value = 0.1
#in orginal model, we cannot regect H0 so there is stationarity, #then we try first difference.in difference Y , we regect H0 so #no stationarity.
par(mfrow=c(2,3))
plot(y, main=expression(y))
plot.new
## function ()
## {
## for (fun in getHook("before.plot.new")) {
## if (is.character(fun))
## fun <- get(fun)
## try(fun())
## }
## .External2(C_plot_new)
## grDevices:::recordPalette()
## for (fun in getHook("plot.new")) {
## if (is.character(fun))
## fun <- get(fun)
## try(fun())
## }
## invisible()
## }
## <bytecode: 0x7f9e1e315780>
## <environment: namespace:graphics>
plot(d.y,main=expression(paste(Delta,"y)")))
#acf and Pacf
acf(coredata(y),type="correlation",lag=40,ylab="",main="ACF")
acf(coredata(y),type="partial",lag=40,ylab="",main="PACF")
acf(coredata(d.y),type="correlation",lag=40,ylab="",main="ACF")
acf(coredata(d.y),type="partial",lag=40,ylab="",main="PACF")
# estimate model -
ARMA1<-arima(y,order=c(1,0,1))
ARMA1
##
## Call:
## arima(x = y, order = c(1, 0, 1))
##
## Coefficients:
## ar1 ma1 intercept
## 0.9906 0.0831 5.5691
## s.e. 0.0045 0.0287 0.7562
##
## sigma^2 estimated as 0.04417: log likelihood = 113.07, aic = -218.14
tsdiag(ARMA1,gof.lag=12)
ARMA2<-arima(y,order=c(2,0,2))
ARMA2
##
## Call:
## arima(x = y, order = c(2, 0, 2))
##
## Coefficients:
## ar1 ar2 ma1 ma2 intercept
## 1.8404 -0.8467 -0.8468 0.2165 5.7271
## s.e. 0.0332 0.0329 0.0465 0.0350 0.3936
##
## sigma^2 estimated as 0.03821: log likelihood = 172.06, aic = -332.12
tsdiag(ARMA2,gof.lag=12)
ARMA3<-arima(y,order=c(3,0,1))
ARMA3
##
## Call:
## arima(x = y, order = c(3, 0, 1))
##
## Coefficients:
## ar1 ar2 ar3 ma1 intercept
## 1.5881 -0.3604 -0.2350 -0.5831 5.7296
## s.e. 0.0604 0.0931 0.0386 0.0555 0.3854
##
## sigma^2 estimated as 0.0383: log likelihood = 171.06, aic = -330.12
tsdiag(ARMA3,gof.lag=12)
ARMA4<-arima(y,order=c(4,0,1))
ARMA4
##
## Call:
## arima(x = y, order = c(4, 0, 1))
##
## Coefficients:
## ar1 ar2 ar3 ar4 ma1 intercept
## 1.5256 -0.3077 -0.1788 -0.0473 -0.5262 5.7201
## s.e. 0.0932 0.1120 0.0697 0.0498 0.0874 0.3884
##
## sigma^2 estimated as 0.03826: log likelihood = 171.5, aic = -329
tsdiag(ARMA4,gof.lag=12)
ARMA5<-arima(y,order=c(3,0,2))
ARMA5
##
## Call:
## arima(x = y, order = c(3, 0, 2))
##
## Coefficients:
## ar1 ar2 ar3 ma1 ma2 intercept
## 2.6080 -2.3318 0.7210 -1.5943 0.8048 5.6760
## s.e. 0.0641 0.1336 0.0706 0.0659 0.0690 0.4967
##
## sigma^2 estimated as 0.03802: log likelihood = 173.85, aic = -333.69
tsdiag(ARMA5,gof.lag=12)
ARMA6<-arima(y,order=c(3,0,3))
ARMA6
##
## Call:
## arima(x = y, order = c(3, 0, 3))
##
## Coefficients:
## ar1 ar2 ar3 ma1 ma2 ma3 intercept
## 0.8706 0.9310 -0.8142 0.1378 -0.6212 0.2290 5.7583
## s.e. 0.0357 0.0232 0.0348 0.0478 0.0453 0.0351 0.3928
##
## sigma^2 estimated as 0.0378: log likelihood = 176.04, aic = -336.09
tsdiag(ARMA6,gof.lag=12)
#based on AIC values and Ljunb-Box tests we choose models 6,4,and2
#forecasting
library("forecast")
## Loading required package: timeDate
## This is forecast 6.2
ARMA2.fcast<-forecast(ARMA2,h=12)
par(mfrow=c(2,2),cex=0.7,mar=c(2,4,3,1))
plot(ARMA2.fcast,xlim=c(2012,2018))
lines(y,12)
ARMA4.fcast<-forecast(ARMA4,h=12)
par(mfrow=c(2,1),cex=0.7,mar=c(2,4,3,1))
plot(ARMA4.fcast,xlim=c(2012,2018))
lines(y,12)
ARMA6.fcast<-forecast(ARMA6,h=12)
par(mfrow=c(2,1),cex=0.7,mar=c(2,4,3,1))
plot(ARMA6.fcast,xlim=c(2012,2018))
lines(y,12)
