ECO 5316 Time Series Econometrics Spring 2016, HW-3,Problem 2

library("zoo")

## Warning: package 'zoo' was built under R version 3.2.3

## 
## Attaching package: 'zoo'

## The following objects are masked from 'package:base':
## 
##     as.Date, as.Date.numeric

library("forecast")

## Warning: package 'forecast' was built under R version 3.2.3

## Loading required package: timeDate

## Warning: package 'timeDate' was built under R version 3.2.3

## This is forecast 6.2

library("tseries")

## Warning: package 'tseries' was built under R version 3.2.3

Problem 2

Consider the monthly time series for Civilian Unemployment Rate, FRED/UNRATE.

Test it for nonstationarity using ADF and KPPS tests.
Follow the Box-Jenkins methodology to build a time series model based on the data until the end of 2014.
After you check the model for adequacy use it to produce and plot a forecast until the end of 2016.

Introduction:

In this problem, data is analysed about monthly Civilian Unemployment Rate from 1948 to 2016.

DATA:

The data for this problem is collectd from Quand.The Civilian Unemployment Rate used in this anlaysis is monthly data from 1948 to 2016.

unrate_d<- read.csv((url("https://research.stlouisfed.org/fred2/data/UNRATE.csv")))

unrate<- ts(unrate_d[,2], start=c(1948,1), frequency = 12)

str(unrate)

##  Time-Series [1:817] from 1948 to 2016: 3.4 3.8 4 3.9 3.5 3.6 3.6 3.9 3.8 3.7 ...

head(unrate)

## [1] 3.4 3.8 4.0 3.9 3.5 3.6

tail(unrate)

## [1] 5.1 5.1 5.0 5.0 5.0 4.9

tsdisplay(unrate, lag.max = 100, xlab=" Years 1948-2016", ylab="ICivilian Unemployment Rate", main="Civilian Unemployment Rate, monthly")

The plot above indicate there is not variance in data set.Lets check for the stationary of the data. This is done conducting Augmented Dickey-Fuller(ADF) and Kwiatkowski-Phillips-Schmidt-Shin (KPSS) test.

adf.test(unrate)

## Warning in adf.test(unrate): p-value smaller than printed p-value

## 
##  Augmented Dickey-Fuller Test
## 
## data:  unrate
## Dickey-Fuller = -4.0267, Lag order = 9, p-value = 0.01
## alternative hypothesis: stationary

kpss.test(unrate, null="Trend")

## Warning in kpss.test(unrate, null = "Trend"): p-value smaller than printed
## p-value

## 
##  KPSS Test for Trend Stationarity
## 
## data:  unrate
## KPSS Trend = 0.49797, Truncation lag parameter = 6, p-value = 0.01

The p-value for ADF and KPSS test is 0.01. This suggests, we reject the null hypothesis,which is non-stationary of the data.

 #split sample - into two parts in order to test the forecasting  adequecy of the model
unrate_all <- unrate
unrate_2014 <- window(unrate_all, end=c(2014,12))

unrate_2015 <- window(unrate_all, start=c(2015,1))

Modeling

Assessing the ACF and PACF plots, let try to identity the ARMA model. ACF is decaying slowely and PACF drops to zero after some lag with some seasonally spike. Counting the lags on PCF plot; PACF drops to zero after 6 lags and 2 or 3 lag seasonal spike around 12 lag period. This suggests AR(6)(2)_12 model.

ar6_3_12 <- arima(unrate_2014, order=c(6,0,0), seasonal=list(order=c(3,0,0),period=12))
tsdisplay(residuals(ar6_3_12))

ar6_3_12

## 
## Call:
## arima(x = unrate_2014, order = c(6, 0, 0), seasonal = list(order = c(3, 0, 0), 
##     period = 12))
## 
## Coefficients:
##          ar1     ar2     ar3      ar4     ar5      ar6     sar1     sar2
##       0.9919  0.1992  -0.031  -0.0525  0.0143  -0.1306  -0.2221  -0.2171
## s.e.  0.0351  0.0499   0.051   0.0503  0.0499   0.0355   0.0365   0.0387
##          sar3  intercept
##       -0.1576     5.7151
## s.e.   0.0387     0.4636
## 
## sigma^2 estimated as 0.03564:  log likelihood = 196.19,  aic = -370.38

tsdiag(ar6_3_12,gof.lag=40)

The above AR(6)(3)_12 model has higher forecasting power for lower lag and the power decrease for higher lag. The forcasting power decrease every years indicating seasonallity.Let’s add MA part in the model.

arma61_31_12 <- arima(unrate_2014, order=c(6,0,1), seasonal=list(order=c(3,0,1),period=12))
tsdisplay(residuals(arma61_31_12))

arma61_31_12

## 
## Call:
## arima(x = unrate_2014, order = c(6, 0, 1), seasonal = list(order = c(3, 0, 1), 
##     period = 12))
## 
## Coefficients:
##          ar1     ar2      ar3      ar4     ar5      ar6      ma1    sar1
##       1.0041  0.1895  -0.0372  -0.0552  0.0232  -0.1309  -0.0166  0.3927
## s.e.  0.2191  0.2224   0.0695   0.0511  0.0544   0.0429   0.2199  0.2120
##          sar2     sar3     sma1  intercept
##       -0.1002  -0.0307  -0.6359     5.7162
## s.e.   0.0603   0.0790   0.2121     0.4699
## 
## sigma^2 estimated as 0.03534:  log likelihood = 199.21,  aic = -372.42

tsdiag(arma61_31_12,gof.lag=70)

The above ARMA(6,1)(3,1)_12 model has higher forecasting power for lower lag and the power decrease for higher lag. The forcasting power decrease every years indicating seasonallity.However, this model is better than previous one.

Forecasting

The forceasting consist of using the second part of the data to understand the forcasting power of the above model.

par(mfrow=c(1,1))
AR.fcast1 <- forecast(arma61_31_12, h=13)
plot(AR.fcast1, xlim=c(2000,2016))
lines(unrate_2015)

As indicated in the modeling section, the forecasting is good for early time period. The forecasting power decreases as the time increases.