HW5
Mai AlSwayan
1st Question
1.Obtain the monthly data for Total Private Residential Construction Spending FRED/PRRESCON. 2.Follow the Box-Jenkins methodology to build a time series model based on the data until the end of 2013. 3.produce and plot a multistep forecast and a sequence of one month ahead forecasts until the end of 2016. 4.Comment on the accuracy of the two forecasts.
Ansewr:
Part 1
Data and data overview
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.numerictprcs <- Quandl("FRED/PRRESCON", type="zoo")
str(tprcs)## 'zooreg' series from Jan 1993 to Feb 2017
## Data: num [1:290] 13430 12141 14007 15762 17390 ...
## Index: Class 'yearmon' num [1:290] 1993 1993 1993 1993 1993 ...
## Frequency: 12summary(tprcs)## Index tprcs
## Min. :1993 Min. :12141
## 1st Qu.:1999 1st Qu.:21414
## Median :2005 Median :27750
## Mean :2005 Mean :29797
## 3rd Qu.:2011 3rd Qu.:35400
## Max. :2017 Max. :60866head(tprcs)## Jan 1993 Feb 1993 Mar 1993 Apr 1993 May 1993 Jun 1993
## 13430 12141 14007 15762 17390 18953tail(tprcs)## Sep 2016 Oct 2016 Nov 2016 Dec 2016 Jan 2017 Feb 2017
## 41204 42111 39135 35397 31847 32225plot(tprcs, xlab="Years", ylab="Million $", main="Total Private Residential Construction Spending ")
Part 2
ltprcs <- log(tprcs)
dltprcs1 <- diff(ltprcs)
dldtprcs <- diff(ltprcs,12)
dlddtprcs <- diff(diff(ltprcs,12))par(mfrow=c(2,3))
plot(tprcs, main=expression(tprcs))
plot(ltprcs, main=expression(log(tprcs)))
plot.new()
plot(dltprcs1, main=expression(paste(Delta, "log(tprcs)")))
plot(dldtprcs, main=expression(paste(Delta[12], "log(tprcs)")))
plot(dlddtprcs, main=expression(paste(Delta, Delta[12], "log(tprcs)")))
library(forecast)maxlag <- 20
par(mfrow=c(2,4), max=c(3,3,4,2))## Warning in par(mfrow = c(2, 4), max = c(3, 3, 4, 2)): "max" is not a
## graphical parameterAcf(ltprcs, type='correlation', lag=maxlag, na.action=na.omit, ylab="", main=expression(paste("ACF for lag(tprcs)")))
Acf(dltprcs1, type='correlation', lag=maxlag, na.action=na.omit, ylab="", main=expression(paste("ACF for",Delta," lag(tprcs)")))
Acf(dldtprcs, type='correlation',lag=maxlag, na.action=na.omit, ylab="", main=expression(paste("ACF for",Delta[12]," lag(tprcs)")))
Acf(dlddtprcs, type='correlation',lag=maxlag, na.action=na.omit, ylab="", main=expression(paste("ACF for",Delta, Delta[12],"lag(tprcs)")))
Acf(ltprcs, type='partial',lag=maxlag, na.action=na.omit, ylab="", main=expression(paste("PACF for lag(tprcs)")))
Acf(dltprcs1, type='partial',lag=maxlag, na.action=na.omit, ylab="", main=expression(paste("PACF for",Delta," lag(tprcs)")))
Acf(dldtprcs, type='partial',lag=maxlag, na.action=na.omit, ylab="", main=expression(paste("PACF for",Delta[12]," lag(tprcs)")))
Acf(dlddtprcs, type='partial',lag=maxlag, na.action=na.omit, ylab="", main=expression(paste("PACF for",Delta,Delta[12]," lag(tprcs)")))
m1 <- Arima(dltprcs1, order = c(1,1,0), seasonal = list(order = c(1,0,0), period = 12))
m1## Series: dltprcs1
## ARIMA(1,1,0)(1,0,0)[12]
##
## Coefficients:
## ar1 sar1
## -0.3996 0.9727
## s.e. 0.0545 0.0090
##
## sigma^2 estimated as 0.000415: log likelihood=696.08
## AIC=-1386.15 AICc=-1386.07 BIC=-1375.16tsdiag(m1, gof.lag=36)
m2 <- Arima(dltprcs1, order = c(1,1,0), seasonal = list(order = c(2,0,0), period = 12))
m2## Series: dltprcs1
## ARIMA(1,1,0)(2,0,0)[12]
##
## Coefficients:
## ar1 sar1 sar2
## -0.4152 0.9164 0.0578
## s.e. 0.0570 0.0669 0.0681
##
## sigma^2 estimated as 0.0004154: log likelihood=696.44
## AIC=-1384.87 AICc=-1384.73 BIC=-1370.22tsdiag(m2, gof.lag=36)
m3 <- Arima(dltprcs1, order = c(2,0,0), seasonal = list(order = c(1,0,0), period = 12))
m3## Series: dltprcs1
## ARIMA(2,0,0)(1,0,0)[12] with non-zero mean
##
## Coefficients:
## ar1 ar2 sar1 mean
## 0.4095 0.2085 0.9715 0.0087
## s.e. 0.0577 0.0577 0.0093 0.0519
##
## sigma^2 estimated as 0.000361: log likelihood=719.88
## AIC=-1429.75 AICc=-1429.54 BIC=-1411.42tsdiag(m3, gof.lag=36)
m4 <- Arima(dltprcs1, order = c(2,0,0), seasonal = list(order = c(1,0,0), period = 12))
m4## Series: dltprcs1
## ARIMA(2,0,0)(1,0,0)[12] with non-zero mean
##
## Coefficients:
## ar1 ar2 sar1 mean
## 0.4095 0.2085 0.9715 0.0087
## s.e. 0.0577 0.0577 0.0093 0.0519
##
## sigma^2 estimated as 0.000361: log likelihood=719.88
## AIC=-1429.75 AICc=-1429.54 BIC=-1411.42tsdiag(m4, gof.lag=36)
m5 <- Arima(dltprcs1, order = c(4,1,0), seasonal = list(order = c(1,0,0), period = 12))
m5## Series: dltprcs1
## ARIMA(4,1,0)(1,0,0)[12]
##
## Coefficients:
## ar1 ar2 ar3 ar4 sar1
## -0.4486 -0.1579 -0.2022 -0.0505 0.9734
## s.e. 0.0590 0.0641 0.0642 0.0595 0.0088
##
## sigma^2 estimated as 0.0004024: log likelihood=701.84
## AIC=-1391.68 AICc=-1391.38 BIC=-1369.7tsdiag(m5, gof.lag=36)
m6 <- Arima(dltprcs1, order = c(2,0,1), seasonal = list(order = c(2,0,0), period = 12))
m6## Series: dltprcs1
## ARIMA(2,0,1)(2,0,0)[12] with non-zero mean
##
## Coefficients:
## ar1 ar2 ma1 sar1 sar2 mean
## -0.3801 0.5466 0.869 0.8735 0.0997 0.0064
## s.e. 0.0638 0.0492 0.057 0.0649 0.0659 0.0439
##
## sigma^2 estimated as 0.0003523: log likelihood=723.89
## AIC=-1433.78 AICc=-1433.38 BIC=-1408.12tsdiag(m6, gof.lag=36)
m7 <- Arima(dltprcs1, order = c(3,0,0), seasonal = list(order = c(1,0,0), period = 12))
m7## Series: dltprcs1
## ARIMA(3,0,0)(1,0,0)[12] with non-zero mean
##
## Coefficients:
## ar1 ar2 ar3 sar1 mean
## 0.4241 0.2402 -0.0734 0.9700 0.0076
## s.e. 0.0587 0.0630 0.0594 0.0098 0.0469
##
## sigma^2 estimated as 0.0003611: log likelihood=720.64
## AIC=-1429.27 AICc=-1428.97 BIC=-1407.27tsdiag(m7, gof.lag=36)
m8 <- Arima(dltprcs1, order = c(9,1,0), seasonal = list(order = c(1,0,0), period = 12))
m8## Series: dltprcs1
## ARIMA(9,1,0)(1,0,0)[12]
##
## Coefficients:
## ar1 ar2 ar3 ar4 ar5 ar6 ar7
## -0.4367 -0.1768 -0.2218 -0.1082 -0.0505 -0.0894 -0.2386
## s.e. 0.0591 0.0647 0.0646 0.0659 0.0667 0.0653 0.0644
## ar8 ar9 sar1
## 0.0146 -0.0221 0.9746
## s.e. 0.0660 0.0615 0.0087
##
## sigma^2 estimated as 0.0003833: log likelihood=710.8
## AIC=-1399.6 AICc=-1398.64 BIC=-1359.31tsdiag(m8, gof.lag=36)
m9 <- Arima(dltprcs1, order = c(9,1,0), seasonal = list(order = c(0,1,0), period = 12))
m9## Series: dltprcs1
## ARIMA(9,1,0)(0,1,0)[12]
##
## Coefficients:
## ar1 ar2 ar3 ar4 ar5 ar6 ar7
## -0.4297 -0.1806 -0.2104 -0.0875 -0.0473 -0.0737 -0.2389
## s.e. 0.0604 0.0659 0.0654 0.0667 0.0676 0.0671 0.0663
## ar8 ar9
## 0.0327 -0.0094
## s.e. 0.0675 0.0620
##
## sigma^2 estimated as 0.0003866: log likelihood=696.97
## AIC=-1373.95 AICc=-1373.12 BIC=-1337.74tsdiag(m9, gof.lag=36)
m10<- Arima(dltprcs1, order = c(9,1,0), seasonal = list(order = c(1,1,0), period = 12))
m10## Series: dltprcs1
## ARIMA(9,1,0)(1,1,0)[12]
##
## Coefficients:
## ar1 ar2 ar3 ar4 ar5 ar6 ar7
## -0.4673 -0.2049 -0.2323 -0.0931 -0.0587 -0.0778 -0.2426
## s.e. 0.0639 0.0683 0.0675 0.0681 0.0691 0.0684 0.0673
## ar8 ar9 sar1
## 0.0177 -0.0247 -0.1219
## s.e. 0.0688 0.0625 0.0681
##
## sigma^2 estimated as 0.0003833: log likelihood=698.56
## AIC=-1375.11 AICc=-1374.11 BIC=-1335.29tsdiag(m10, gof.lag=36)
m8 is the best model so far to use for our forecast.
Part 3
end.month <- 2013 + 11/12
tprcs1 <- window( tprcs, end = end.month)
tprcs2 <- window( tprcs, start = end.month + 1/12)Multistep Forecast
multiple<- forecast(m8, length(tprcs2))
multiple## Point Forecast Lo 80 Hi 80 Lo 95
## Mar 2017 0.184837190 0.159745817 0.2099285631 0.146463246
## Apr 2017 0.036742960 0.007944035 0.0655418839 -0.007301195
## May 2017 0.043945776 0.011710629 0.0761809221 -0.005353627
## Jun 2017 0.068977890 0.035012084 0.1029436967 0.017031672
## Jul 2017 0.039489631 0.003463913 0.0755153491 -0.015606951
## Aug 2017 0.003965373 -0.033989917 0.0419206622 -0.054082233
## Sep 2017 -0.032577277 -0.072053684 0.0068991289 -0.092951231
## Oct 2017 0.027250963 -0.012746705 0.0672486310 -0.033920192
## Nov 2017 -0.067453383 -0.109091098 -0.0258156676 -0.131132773
## Dec 2017 -0.088062481 -0.130845615 -0.0452793479 -0.153493638
## Jan 2018 -0.098251604 -0.142636810 -0.0538663975 -0.166132918
## Feb 2018 0.016754778 -0.028806304 0.0623158602 -0.052924884
## Mar 2018 0.184407788 0.126756594 0.2420589828 0.096237895
## Apr 2018 0.043660256 -0.019089823 0.1064103349 -0.052307709
## May 2018 0.047916755 -0.019910332 0.1157438424 -0.055815823
## Jun 2018 0.073655691 0.002771702 0.1445396804 -0.034752015
## Jul 2018 0.043013283 -0.031443876 0.1174704424 -0.070859115
## Aug 2018 0.010543983 -0.067074909 0.0881628760 -0.108163869
## Sep 2018 -0.026088231 -0.106677972 0.0545015103 -0.149339604
## Oct 2018 0.032735707 -0.049599581 0.1150709941 -0.093185249
## Nov 2018 -0.060583633 -0.145821324 0.0246540577 -0.190943432
## Dec 2018 -0.079613390 -0.167211475 0.0079846948 -0.213583100
## Jan 2019 -0.090207488 -0.180651543 0.0002365665 -0.228529733
## Feb 2019 0.022546590 -0.070228774 0.1153219536 -0.119341085
## Mar 2019 0.185140231 0.081868003 0.2884124584 0.027198988
## Apr 2019 0.048552825 -0.060412554 0.1575182037 -0.118095341
## May 2019 0.052363378 -0.062387263 0.1671140181 -0.123132583
## Jun 2019 0.077810941 -0.040998466 0.1966203485 -0.103892368
## Jul 2019 0.047552066 -0.075763172 0.1708673046 -0.141042317
## Aug 2019 0.016203636 -0.111213687 0.1436209593 -0.178664344
## Sep 2019 -0.019752267 -0.151115022 0.1116104876 -0.220654264
## Oct 2019 0.037865867 -0.096274083 0.1720058177 -0.167283484
## Nov 2019 -0.053343419 -0.191323307 0.0846364688 -0.264365447
## Dec 2019 -0.071709437 -0.212999709 0.0695808342 -0.287794261
## Jan 2020 -0.082188134 -0.227223027 0.0628467584 -0.303999861
## Feb 2020 0.027858470 -0.120429328 0.1761462693 -0.198928147
## Mar 2020 0.186196443 0.028033961 0.3443589245 -0.055692200
## Apr 2020 0.053180154 -0.111157797 0.2175181053 -0.198153054
## Hi 95
## Mar 2017 0.223211134
## Apr 2017 0.080787114
## May 2017 0.093245178
## Jun 2017 0.120924108
## Jul 2017 0.094586212
## Aug 2017 0.062012979
## Sep 2017 0.027796676
## Oct 2017 0.088422118
## Nov 2017 -0.003773993
## Dec 2017 -0.022631325
## Jan 2018 -0.030370289
## Feb 2018 0.086434440
## Mar 2018 0.272577682
## Apr 2018 0.139628220
## May 2018 0.151649334
## Jun 2018 0.182063398
## Jul 2018 0.156885681
## Aug 2018 0.129251836
## Sep 2018 0.097163142
## Oct 2018 0.158656662
## Nov 2018 0.069776166
## Dec 2018 0.054356320
## Jan 2019 0.048114757
## Feb 2019 0.164434265
## Mar 2019 0.343081473
## Apr 2019 0.215200991
## May 2019 0.227859338
## Jun 2019 0.259514250
## Jul 2019 0.236146449
## Aug 2019 0.211071616
## Sep 2019 0.181149730
## Oct 2019 0.243015218
## Nov 2019 0.157678609
## Dec 2019 0.144375386
## Jan 2020 0.139623593
## Feb 2020 0.254645088
## Mar 2020 0.428085085
## Apr 2020 0.304513362plot(multiple, type="o", pch=16, xlim=c(2013,2016), ylim=c(-0.03,0.03),
main="Multistep Forecasts")
lines(multiple$mean, type="p", pch=16, lty="dashed", col="blue")
lines(dltprcs1, type="o", pch=16, lty="dashed")
we can see the both the actual model and the forcasted model has similar results.