#install.packages('Quandl')
library(Quandl)
Quandl.api_key("XzdSwkDsE98Mxj3ixQzG")
rGDP <- Quandl("FRED/GDPC1", type="zoo" )
str(rGDP)

#The plot for Real GDP shows the upward sloping trend
plot(rGDP, xlab="Year", ylab="Real GDP", main="The Graph of Real GDP")

#Construct the time series with log changes in Real GDP ( ) and denote the variable name for by dlrGDP
dlrGDP <- diff(log(rGDP))

#The plot shows that the volatiliy of real GDP is decreasing over time and we need to check stationariy later
plot(dlrGDP, xlab="Year", ylab="log difference in Real GDP", main="Log Differences in Real
GDP")


#Split the sample into two parts: up to 2008.Q4 vs from 2009.Q1
first.q <- 1947.25
last.q <-2008.75
dlrGDP.p1 <- window(dlrGDP, end=last.q)
dlrGDP.p2 <- window(dlrGDP, start=last.q+0.25)
# We can also split the sample period using the following method.
# dlrGDP.p1 <- window(dlrGDP, end="2008 Q4")
# dlrGDP.p2 <- window(dlrGDP, start="2009 Q1")

#install.packages('forecast')
library(forecast)
#We use ACF and PACF to identify suitable time series model.If ACF dies out slowly and PACF drops to zero
#suddenly after lag p, then it is AR(p) model. If ACF drops to zero immediately after lag q and PACF dies out
#slowly, then it is MA(q) model
Acf(dlrGDP.p1, type="correlation", lag=48, main="ACF for Real GDP")
Acf(dlrGDP.p1, type="partial", lag=48, main="PACF for Real GDP")

#non seasonal data
a1 <- auto.arima(dlrGDP.p1, seasonal = FALSE, stationary = TRUE, stepwise = FALSE, ic ="aicc")
a1

#seasonal data
a2 <- auto.arima(dlrGDP.p1, seasonal = TRUE, stationary = TRUE, stepwise = FALSE, ic ="aicc")
a2


library(forecast)
arma22 <- Arima(dlrGDP.p1, order=c(2, 0, 2))
#plot.Arima(arma22)

#Diagnose residuals
tsdiag(arma22, gof.lag=36)

#The multistep forecasts using ARMA(2,2)
arma22 <- Arima(dlrGDP.p1, order=c(2, 0, 2))
arma22.f <- forecast(arma22, length(dlrGDP.p2))
plot(arma22.f, type="o", pch=16, xlim=c(2005, 2016), ylim=c(-0.03, 0.04), main="ARMA(2,2) Mod
el Multistep Forecasts")
lines(arma22.f$mean, type="p", pch=16, lty="dashed", col="blue")
lines(dlrGDP.p2, type="o", pch=16, lty="dotted", col = "red")

#The one-step rolling forecats
rol.f <-zoo()
for(i in 1: length(dlrGDP.p2))
 {
   y <- window(dlrGDP, end = last.q+(i-1)/4)
   rol.updt <- arima(y, order=c(2,0,2))
   rol.f <- c(rol.f, forecast(rol.updt, 1)$mean)
 } 
rol.f <-as.ts(rol.f)
plot(rol.f, type="o", pch=15, xlim=c(2005, 2016), ylim=c(-0.03, 0.04), main="One Step Rolling
Scheme Forecasts")
lines(rol.f, type="p", pch=15, lty="solid", col="red")
lines(dlrGDP, type="o", pch=16, lty="dotted")
lines(dlrGDP.p1, type="o", pch=16, lty="solid")

#Comparison: Multistep forecast using ARMA(2,2) vs. One step rolling scheme forecast
plot(arma22.f, type="o", pch=16, xlim=c(2005, 2016), ylim=c(-0.03, 0.03), main="Multistep For
ecasts,ARMA(2,2) vs. One step Rolling Scheme Forecast")
lines(arma22.f$mean, type="p", pch=16, lty="dashed", col="blue")
lines(rol.f, type="b", pch=15, lty="solid", lwd=2, col="red")
lines(dlrGDP, type="o", pch=16, lty="dotted")

# Comparison: Multistep forecasts using AR(1) vs. Multistep forecasts using ARMA(2,2)
ar1 <- Arima(dlrGDP.p1, order=c(1, 0, 0))
ar1.f <- forecast(ar1, length(dlrGDP.p2))
plot(ar1.f, type="o", pch=16, xlim=c(2005, 2016), ylim=c(-0.03, 0.04), main="Multistep Foreca
sts,ARMA(2,2) vs. Multistep Forecasts,AR(1)")
lines(ar1.f$mean, type="p", pch=15, lty="dashed", col="green")
lines(arma22.f$mean, type="p", pch=16, lty="dashed", col="blue")
lines(dlrGDP, type="o", pch=16, lty="dotted")

accuracy(arma22.f$mean, dlrGDP.p2)

accuracy(rol.f, dlrGDP.p2)