Importing packages and reading in CSVs


library(forecast)
library(tidyverse)
library(base)

PCE <- read.csv("C:/Users/qures/Downloads/PCEPII.csv",header=TRUE)

PCE and Inflation Series-Stationarity, First Difference, and Structural Break

January 1959 to December 2019

PCE <-  subset(PCE, DATE>="1959-01-01" & DATE <="2019-12-01")
PCE$DATE <- as.Date(DATE,format = "%Y-%m-%d")
PCE$LogPCEPI <- log(PCEPI)
PCE$LoglPCEPI <- lag(PCE$LogPCEPI)
PCE$infl <- 1200*(PCE$LogPCEPI-PCE$LoglPCEPI)

plot(PCE$DATE,PCE$PCEPI,type="l",main="PCEPI Time Series",xlab = "Date", ylab = "PCEPI")
mtext("Jan 1959-Dec 2019", side = 3, line = 2.5)

Based on visual inspection, the PCEPI series is not stationary as its plot shows an exponential trend and the mean and variance seem to depend on the previous value.

plot(PCE$DATE,PCE$infl,type = "l",main="Annualized Inflation Rate",xlab = "Date",ylab="Inflation")
mtext("Jan 1959-Dec 2019", side = 3, line = 2.5)

par(mfrow=c(1,2))
ACF <- acf(PCE$infl, lag.max = NULL, type = c("correlation"),plot = FALSE, na.action = na.pass,demean = TRUE)
plot(ACF,main="Jan 1959 - Dec 2019")
mtext("Inflation Series", side = 3, line = 3)
PACF <- acf(PCE$infl, lag.max = NULL, type = c("partial"),plot = FALSE, na.action = na.pass,demean = TRUE)
plot(PACF,main="Jan 1959 - Dec 2019")
mtext("Inflation Series", side = 3, line = 3)

After inspection of the ACF, it seems that the annualized inflation series is non-stationary as the ACF is decreasing for the lags and it is a positive number. It does not seem to be suitable for a finite ARMA model because the plot shows persistence and the ACF does not converge towards 0, which shows that it is not suitable for a finite ARMA model.

deltinfl <- diff(PCE$infl)

par(mfrow=c(1,1))
plot(deltinfl,type="l",main="Change in Inflation")
mtext("Jan 1959-Dec 2019", side = 3, line = 2.5)

par(mfrow=c(1,2))
ACF_diff <- acf(deltinfl,lag.max = NULL, type = c("correlation"),plot = FALSE, na.action = na.pass,demean = TRUE)
plot(ACF_diff,main="Jan 1959- Dec 2019")
mtext("Change in Inflation Series", side = 3, line = 3)
PACF_diff <- acf(deltinfl,lag.max = NULL, type = c("partial"),plot = FALSE, na.action = na.pass,demean = TRUE)
plot(PACF_diff,main="Jan 1959- Dec 2019")
mtext("Change in Inflation Series", side = 3, line = 3)


The inflation change growth series appears to be stationary after visual inspection and it is suitable for an MA(2) model.

Pre-1984 (January 1959 to December 1984) 

PC84 <-  subset(PCE, DATE>="1959-01-01" & DATE <="1984-12-01")


plot(PC84$DATE,PC84$PCEPI,type = "l",main="PCEPI Time Series",xlab="Date",ylab="PCEPI")
mtext("Jan 1959-Dec 1984", side = 3, line = 2.5)

The PCEPI series the sub-sample of January 1959 to December 1984 appears to be nonstationary and it also has an exponential trend.

plot(PC84$DATE,PC84$infl,type = "l",main="Annualized Inflation Rate",xlab="Date",ylab="Inflation")
mtext("Jan 1959-Dec 1984", side = 3, line = 2.5)

par(mfrow=c(1,2))
ACF84 <- acf(PC84$infl, lag.max = NULL, type = c("correlation"),plot = FALSE, na.action = na.pass,demean = TRUE)
plot(ACF84,main="Jan 1959 - Dec 1984")
mtext("Inflation Series", side = 3, line = 3)

PACF84 <- acf(PC84$infl, lag.max = NULL, type = c("partial"),plot = FALSE, na.action = na.pass,demean = TRUE)
plot(PACF84,main="Jan 1959 - Dec 1984")
mtext("Inflation Series", side = 3, line = 3)

The inflation time series for the sub-sample of January 1959 to December 1984 is not stationary and it is also not suitable for an ARMA model.

diffinfl84 <- diff(PC84$infl)
par(mfrow=c(1,1))
plot(diffinfl84, type = "l", main = "Change in inflation for 1959-1984 ")

par(mfrow=c(1,2))
ACF84_diff <- acf(diffinfl84,lag.max = NULL, type = c("correlation"),plot = FALSE, na.action = na.pass,demean = TRUE)
plot(ACF84_diff,main="Jan 1959 - Dec 1984")
mtext("Change in Inflation Series", side = 3, line = 3)

PACF84_diff <- acf(diffinfl84,lag.max = NULL, type = c("partial"),plot = FALSE, na.action = na.pass,demean = TRUE)
plot(PACF84_diff,main="Jan 1959 - Dec 1984")
mtext("Change in Inflation Series", side = 3, line = 3)

The inflation growth series appears to be stationary, and it does look suitable for an MA(1) model.



cat("The average annualized monthly inflation over the 1959-1984 subsample is ", mean(PC84$infl, na.rm = TRUE),"\n\n")
## The average annualized monthly inflation over the 1959-1984 subsample is  4.531957



Post-1984 (January 1985 to December 2019) 

PC85 <-  subset(PCE, DATE>="1985-01-01" & DATE <="2019-12-01")
plot(PC85$DATE,PC85$PCEPI,type = "l",main="PCEPI Time Series",xlab="Date",ylab="PCEPI")
mtext("Jan 1985-Dec 2019", side = 3, line = 2.5)

par(mfrow=c(1,1))
plot(PC85$DATE,PC85$infl,type = "l",main="Annualized Inflation Rate",xlab="Date",ylab="Inflation")
mtext("Jan 1985-Dec 2019", side = 3, line = 2.5)

The PCEPI series the sub-sample of January 1985 to December 2019 appears to be nonstationary and it also has an exponential trend

par(mfrow=c(1,2))
ACF85 <- acf(PC85$infl, lag.max = NULL, type = c("correlation"),plot = FALSE, na.action = na.pass,demean = TRUE)
plot(ACF85,main="Jan 1985- Dec 2019")
mtext("Inflation Series", side = 3, line = 3)
PACF85 <- acf(PC85$infl, lag.max = NULL, type = c("partial"),plot = FALSE, na.action = na.pass,demean = TRUE)
plot(PACF85,main="Jan 1985- Dec 2019")
mtext("Inflation Series", side = 3, line = 3)

The inflation time series for the sub-sample of January 1985 to December 2019 is not stationary but does seem suitable for an MA(4) model.

diffinfl85 <- diff(PC85$infl)
par(mfrow=c(1,1))
plot(diffinfl85, type = "l", main = "Change in inflation for 1985-2019 ")


The inflation growth series appears to be stationary, and it does look suitable for an MA(2) model.

par(mfrow=c(1,2))
ACF85_diff <- acf(diffinfl85,lag.max = NULL, type = c("correlation"),plot = FALSE, na.action = na.pass,demean = TRUE)
plot(ACF85_diff,main="Jan 1985- Dec 2019")
mtext("Change in Inflation Series", side = 3, line = 3)


PACF85_diff <- acf(diffinfl85,lag.max = NULL, type = c("partial"),plot = FALSE, na.action = na.pass,demean = TRUE)
plot(PACF85_diff,main="Jan 1985- Dec 2019")
mtext("Change in Inflation Series", side = 3, line = 3)

cat("The average annualzed monthly inflation over the 1985-2019 subsample is ", mean(PC85$infl, na.rm = TRUE),"\n\n")
## The average annualzed monthly inflation over the 1985-2019 subsample is  2.165385
print(t.test(PC84$infl,PC85$infl))
## 
##  Welch Two Sample t-test
## 
## data:  PC84$infl and PC85$infl
## t = 11.037, df = 521.5, p-value < 2.2e-16
## alternative hypothesis: true difference in means is not equal to 0
## 95 percent confidence interval:
##  1.945327 2.787817
## sample estimates:
## mean of x mean of y 
##  4.531957  2.165385
cat("There does seem to be a difference between the means.  \n\n")
## There does seem to be a difference between the means.


For the pre-1984 period, an MA model is not suitable for the inflation series but it is for the change in the inflation series, with an order of 1, i.e., MA(1). But for the post-1984 period, an MA(4) model is suitable for the inflation series and a MA(2) for the change in inflation series. As we can see here, the orders of the models changed for the periods.
For the pre-1984 period, an AR(5) or possibly an AR(6) model is suitable for the inflation series and an AR(4) or possibly an AR(5) is suitable for the growth in inflation series. For the post-1984 period, an AR(1) is possibly suitable for the inflation series an AR(7) or possibly an AR(8) is suitable for the growth in inflation series. There is a change in orders in the two periods.


AR(1) Forecast for Post-1984 Inflation Change Series

diffinfl85_lag <- lag(diffinfl85, n= 1, default = NA)
m1 <- lm(diffinfl85 ~ diffinfl85_lag)
coef(m1)
##    (Intercept) diffinfl85_lag 
##    -0.01054197    -0.24870163
m_pred <- -0.01054 - 0.2470*1.42981974
m_pred
## [1] -0.3637055
cat("The hand-coded prediction is",m_pred,". \n\n")
## The hand-coded prediction is -0.3637055 .
m2 <- predict(ar(diffinfl85,demean=TRUE,order.max=1),n.ahead=1)
m2[["pred"]]
## Time Series:
## Start = 420 
## End = 420 
## Frequency = 1 
## [1] -0.366011
cat("The predict() prediction is", m2[["pred"]],". \n\n")
## The predict() prediction is -0.366011 .
m3 <- arima(diffinfl85,order=c(1,0,0))
forec <- forecast(m3,h=1)
forec[["mean"]]
## Time Series:
## Start = 420 
## End = 420 
## Frequency = 1 
## [1] -0.3660654
cat("The forecast() prediction is",forec[["mean"]],". \n\n")
## The forecast() prediction is -0.3660654 .

Looking at these results, we can see that the forecasted values obtained in the three ways are close to each other.

#Method 1
msem1<- sum(m1$residuals^2)/(418-2)
cat("The MSE using Method 1 is",msem1,". \n\n")
## The MSE using Method 1 is 5.383023 .
#Method 2
msem2<- (421/417)* (sum(resid(m1)^2)/419)
cat("The MSE using Method 2 is",msem2,". \n\n")
## The MSE using Method 2 is 5.395747 .
#Method 3
error <- numeric(84)
forecast_meth3 <- numeric(84)
forecasts = list(length=84)
for (i in 1:84){
  lm2 <- lm(diffinfl85[2:(334+i)]~diffinfl85[1:(333+i)])
  forecast_meth3[i] <- lm2$coef[1]+lm2$coef[2]*diffinfl85[335+i]
  error[i] <- diffinfl85[335+i]-forecast_meth3
}
P20<-(1/84)*sum(error^2)

cat("The MSE using Method 3 is",P20,". \n\n")
## The MSE using Method 3 is 4.289948 .
#predict()
msepred<-m2[["se"]]^2
cat("The MSE using pred() is",msepred,". \n\n")
## The MSE using pred() is 5.370871 .
#forecast()
m3fc<-forec[["model"]][["sigma2"]]
cat("The MSE using forecast() is",m3fc,". \n\n")
## The MSE using forecast() is 5.344903 .

From these results, we can see that the closest to the predict() and forecast() seems to be Method 1, with Method 2 not being too far off from their value. Method 3 seems to have given a value that was quite different from the rest, with a lower value.

Creating lag variables and splitting the data



ts_pcepi<-ts(PCE$PCEPI[2:732])
ts_lagpcepi<-ts(PCE$PCEPI[1:731])
ts_inflation<-1200*(log(ts_pcepi)-log(ts_lagpcepi))
ts_inflation_change<-diff(ts_inflation)
daterange <- ts_inflation_change[418:730]
daterange2 <-ts_inflation_change[311:730]



p19 <- ts(PCE$PCEPI[314:732])
lagp19<- ts(PCE$PCEPI[313:731])
tsinfl <- 1200*(log(p19)-log(lagp19))
fdinfl<- diff(tsinfl)



MSFEpoos based on different out-of-sample length

#P=10%

error1<- numeric(42)
forecast_meth31 <- numeric(42)


for (i in 1:42){
  lm1 <- lm(fdinfl[2:(375+i)]~fdinfl[1:(374+i)])
  forecast_meth31[i] <- lm1$coef[1]+lm1$coef[2]*fdinfl[376+i]
  error1[i] <- fdinfl[376+i]-forecast_meth31
}
P10<-(1/42)*sum(error1^2)
P10
## [1] 3.563745
#P=20%

error2 <- numeric(84)
forecast_meth32 <- numeric(84)

for (i in 1:84){
  lm2 <- lm(fdinfl[2:(333+i)]~fdinfl[1:(332+i)])
  forecast_meth32[i] <- lm2$coef[1]+lm2$coef[2]*fdinfl[334+i]
  error2[i] <- fdinfl[334+i]-forecast_meth32
}
P20<-(1/84)*sum(error2^2)
P20
## [1] 4.286859
#P=30%
error3<- numeric(125)
forecast_meth33 <- numeric(125)

for (i in 1:125){
  lm3 <- lm(fdinfl[2:(292+i)]~fdinfl[1:(291+i)])
  forecast_meth33[i] <- lm3$coef[1]+lm3$coef[2]*fdinfl[293+i]
  error3[i] <- fdinfl[293+i]-forecast_meth33
}
P30<-(1/125)*sum(error3^2)
P30
## [1] 4.437319
#40%
error4<- numeric(167)
forecast_meth34 <- numeric(167)

for (i in 1:167){
  lm4 <- lm(fdinfl[2:(250+i)]~fdinfl[1:(249+i)])
  forecast_meth34[i] <- lm4$coef[1]+lm4$coef[2]*fdinfl[251+i]
  error4[i] <- fdinfl[251+i]-forecast_meth34
}
P40<-(1/167)*sum(error4^2)
P40
## [1] 8.379489
#50%
error5<- numeric(209)
forecast_meth35 <- numeric(209)

for (i in 1:209){
  lm5 <- lm(fdinfl[2:(208+i)]~fdinfl[1:(207+i)])
  forecast_meth35[i] <- lm5$coef[1]+lm5$coef[2]*fdinfl[209+i]
  error5[i] <- fdinfl[209+i]-forecast_meth35
}
P50<-(1/209)*sum(error5^2)
P50
## [1] 6.647351
#60%
error6<- numeric(251)
forecast_meth36<- numeric(251)

for (i in 1:251){
  lm6 <- lm(fdinfl[2:(166+i)]~fdinfl[1:(165+i)])
  forecast_meth36[i] <- lm2$coef[1]+lm2$coef[2]*fdinfl[167+i]
  error6[i] <- fdinfl[167+i]-forecast_meth36
}
P60<-(1/251)*sum(error6^2)
P60
## [1] 7.33955
#70%
error7<- numeric(293)
forecast_meth37 <- numeric(293)

for (i in 1:293){
  lm7 <- lm(fdinfl[2:(124+i)]~fdinfl[1:(123+i)])
  forecast_meth37[i] <- lm7$coef[1]+lm7$coef[2]*fdinfl[125+i]
  error7[i] <- fdinfl[125+i]-forecast_meth37
}
P70<-(1/293)*sum(error7^2)
P70
## [1] 6.533051
#80%
error8<- numeric(334)
forecast_meth38 <- numeric(334)

for (i in 1:334){
  lm8 <- lm(fdinfl[2:(83+i)]~fdinfl[1:(82+i)])
  forecast_meth38[i] <- lm8$coef[1]+lm8$coef[2]*fdinfl[84+i]
  error8[i] <- fdinfl[84+i]-forecast_meth38
}
P80<-(1/334)*sum(error8^2)
P80
## [1] 5.858293

From the results, we see that the MSFEpoos value increases with an increase in P, but then decreases to 5.86 for when P=80%, which is the closest to the value of the MSFEfpe

Back-testing AR(1) model


#AR(1) Model

realfcst<-numeric()
msfe<-numeric()
lowerCI<-numeric()
upperCI<-numeric()


for (i in 1:312){
  lm1<-lm(ts_inflation_change[299:417+i]~ts_inflation_change[298:416+i])
  realfcst[i]<-lm1$coef[1] + lm1$coef[2]*ts_inflation_change[417+i]
  msfe[i]<-((119+i+1+1)/(119+i-1-1))*(sum(lm1$residuals^2)/(119+i))
  lowerCI[i]<-realfcst[i]-1.96*sqrt(msfe[i])
  upperCI[i]<-realfcst[i]+1.96*sqrt(msfe[i])
}

autoplot(ts(realfcst),main="Back-testing on AR(1)",color='blue',ylab='Forecast values')+autolayer(ts(lowerCI),color='red')+autolayer(ts(upperCI),color='red')+
  autolayer(ts(daterange),color='black')

CI<-logical()
for (i in 1:312){
  if (ts_inflation_change[417+i] > lowerCI[i] & ts_inflation_change[417+i] < upperCI[i]){
    CI[i]=TRUE
  }else{
    CI[i]=FALSE
  }
}
sum(CI)
## [1] 228
cat("The realizations were in the confidence intervals,",(sum(CI)/312)*100,"% of the time. 
The mean sum of squares was",(1/312)*sum((abs(realfcst-ts_inflation_change[419:730]))^2),"")
## The realizations were in the confidence intervals, 73.07692 % of the time. 
## The mean sum of squares was 5.944084



Back-testing AR(2) model



#AR(2) Model

realfcst2<-numeric()
msfe2<-numeric()
lowerCI2<-numeric()
upperCI2<-numeric()

for (i in 1:312){
  lm2<-lm(ts_inflation_change[300:417+i]~ts_inflation_change[299:416+i]+ts_inflation_change[298:415+i])
  realfcst2[i]<-lm2$coef[1] + lm2$coef[2]*ts_inflation_change[417+i] + lm2$coef[3]*ts_inflation_change[416+i]
  msfe2[i]<-((119+i+2+1)/(119+i-2-1))*(sum(lm2$residuals^2)/(119+i))
  lowerCI2[i]<-realfcst2[i]-1.96*sqrt(msfe2[i])
  upperCI2[i]<-realfcst2[i]+1.96*sqrt(msfe2[i])
}

autoplot(ts(realfcst2),main="Back-testing on AR(2)",color='blue')+autolayer(ts(lowerCI2),color='red')+autolayer(ts(upperCI2),color='red')+
  autolayer(ts(daterange),color='black')

CI2<-logical()

for (i in 1:312){
  if (ts_inflation_change[417+i] > lowerCI2[i] & ts_inflation_change[417+i] < upperCI2[i])
  {
    CI2[i]=TRUE
  }
  else
  {
    CI2[i]=FALSE
  }
}
sum(CI2)
## [1] 215
cat("The realizations were in the confidence intervals,",(sum(CI2)/312)*100,"% of the time. 
The mean sum of squares was",(1/312)*sum((abs(realfcst2-ts_inflation_change[419:730]))^2),"")
## The realizations were in the confidence intervals, 68.91026 % of the time. 
## The mean sum of squares was 5.620076



Back-testing ARMA(1,2) model

#ARMA(1,2) model

realfcst3<-numeric()
msfe3<-numeric()
lowerCI3<-numeric()
upperCI3<-numeric()

for (i in 1:312){
  fsct<-forecast(arima(ts_inflation_change[300:417+i], order = c(1,0,2)),h=1)
  realfcst3[i]<-fsct[["mean"]][1]
  msfe3[i]<-fsct[["model"]][["sigma2"]]
  lowerCI3[i]<-realfcst3[i]-1.96*sqrt(msfe3[i])
  upperCI3[i]<-realfcst3[i]+1.96*sqrt(msfe3[i])
}

autoplot(ts(realfcst3),main="Back-testing on ARMA(1,2)",,color='blue',ylab='Forecast values')+autolayer(ts(lowerCI3),color='red')+
  autolayer(ts(upperCI3),color='red')+autolayer(ts(daterange),color='black')

CI3<-logical()
for (i in 1:312)
{
  if (ts_inflation_change[417+i] > lowerCI3[i] & ts_inflation_change[417+i] < upperCI3[i]){
    CI3[i]=TRUE
  }
  else
  {
    CI3[i]=FALSE
  }
}
sum(CI3)
## [1] 241
cat("The realizations were in the confidence intervals,",(sum(CI3)/312)*100,"% of the time. 
The mean sum of squares was",(1/312)*sum((abs(realfcst3-ts_inflation_change[419:730]))^2),"")
## The realizations were in the confidence intervals, 77.24359 % of the time. 
## The mean sum of squares was 4.582007



Model Selection


AR1 <- arima(ts_inflation_change,order=c(1,0,0)) 
AR2 <- arima(ts_inflation_change,order=c(2,0,0))
ARMA <- arima(ts_inflation_change,order=c(1,0,2))
MA3 <- arima(ts_inflation_change,order=c(0,0,3)) 

AR1$aic
## [1] 3204.003
AR2$aic
## [1] 3160.906
ARMA$aic
## [1] 3060.71
MA3$aic
## [1] 3059.972
auto<-auto.arima(ts_inflation_change)
auto$aic
## [1] 3057.972

Single Forecast Over Multiple Horizons


model <- arima(ts_inflation_change[311:730],order=c(1,0,0))
focst <- forecast(model,h=3)
summary(focst)
## 
## Forecast method: ARIMA(1,0,0) with non-zero mean
## 
## Model Information:
## 
## Call:
## arima(x = ts_inflation_change[311:730], order = c(1, 0, 0))
## 
## Coefficients:
##           ar1  intercept
##       -0.2486    -0.0043
## s.e.   0.0473     0.0904
## 
## sigma^2 estimated as 5.347:  log likelihood = -948.07,  aic = 1902.13
## 
## Error measures:
##                       ME     RMSE     MAE      MPE     MAPE      MASE       ACF1
## Training set 0.001354483 2.312405 1.72315 81.28657 133.9653 0.5960368 -0.0710918
## 
## Forecasts:
##     Point Forecast     Lo 80    Hi 80     Lo 95    Hi 95
## 421    -0.36074807 -3.324214 2.602718 -4.892979 4.171483
## 422     0.08435434 -2.969300 3.138009 -4.585807 4.754516
## 423    -0.02628998 -3.085430 3.032850 -4.704841 4.652261
model2 <- arima(ts_inflation_change[311:730],order=c(2,0,0))
focst2 <- forecast(model2,h=3)
summary(focst2)
## 
## Forecast method: ARIMA(2,0,0) with non-zero mean
## 
## Model Information:
## 
## Call:
## arima(x = ts_inflation_change[311:730], order = c(2, 0, 0))
## 
## Coefficients:
##           ar1      ar2  intercept
##       -0.3196  -0.2865    -0.0050
## s.e.   0.0468   0.0467     0.0674
## 
## sigma^2 estimated as 4.906:  log likelihood = -930.07,  aic = 1868.15
## 
## Error measures:
##                      ME     RMSE      MAE      MPE     MAPE      MASE       ACF1
## Training set 0.00254465 2.214983 1.615356 53.48843 133.8038 0.5587509 -0.0643723
## 
## Forecasts:
##     Point Forecast     Lo 80    Hi 80     Lo 95    Hi 95
## 421     0.08070127 -2.757913 2.919316 -4.260585 4.421987
## 422    -0.44347184 -3.423576 2.536632 -5.001148 4.114204
## 423     0.11063824 -2.915064 3.136341 -4.516774 4.738051
model3 <- arima(ts_inflation_change[311:730],order=c(1,0,2))
focst3 <- forecast(model3,h=3)
summary(focst3)
## 
## Forecast method: ARIMA(1,0,2) with non-zero mean
## 
## Model Information:
## 
## Call:
## arima(x = ts_inflation_change[311:730], order = c(1, 0, 2))
## 
## Coefficients:
##          ar1      ma1      ma2  intercept
##       0.2398  -0.8080  -0.1920    -0.0047
## s.e.  0.1130   0.1141   0.1137     0.0012
## 
## sigma^2 estimated as 3.96:  log likelihood = -887.66,  aic = 1785.31
## 
## Error measures:
##                     ME     RMSE      MAE      MPE     MAPE      MASE         ACF1
## Training set -0.040883 1.990098 1.441103 44.32186 146.8271 0.4984767 -0.001788268
## 
## Forecasts:
##     Point Forecast     Lo 80    Hi 80     Lo 95    Hi 95
## 421     -0.6505800 -3.204009 1.902849 -4.555712 3.254552
## 422     -0.4385101 -3.372401 2.495381 -4.925510 4.048489
## 423     -0.1087294 -3.159273 2.941814 -4.774133 4.556674
autoplot(ts(focst[["mean"]]),main="3 month forecasts")+autolayer(ts(focst2[["mean"]]))+autolayer(ts(focst3[["mean"]]))

#####

model4 <- arima(ts_inflation_change[311:730],order=c(1,0,0))
focst4 <- forecast(model4,h=6)
summary(focst4)
## 
## Forecast method: ARIMA(1,0,0) with non-zero mean
## 
## Model Information:
## 
## Call:
## arima(x = ts_inflation_change[311:730], order = c(1, 0, 0))
## 
## Coefficients:
##           ar1  intercept
##       -0.2486    -0.0043
## s.e.   0.0473     0.0904
## 
## sigma^2 estimated as 5.347:  log likelihood = -948.07,  aic = 1902.13
## 
## Error measures:
##                       ME     RMSE     MAE      MPE     MAPE      MASE       ACF1
## Training set 0.001354483 2.312405 1.72315 81.28657 133.9653 0.5960368 -0.0710918
## 
## Forecasts:
##     Point Forecast     Lo 80    Hi 80     Lo 95    Hi 95
## 421   -0.360748069 -3.324214 2.602718 -4.892979 4.171483
## 422    0.084354337 -2.969300 3.138009 -4.585807 4.754516
## 423   -0.026289979 -3.085430 3.032850 -4.704841 4.652261
## 424    0.001214174 -3.058265 3.060693 -4.677855 4.680283
## 425   -0.005622855 -3.065123 3.053877 -4.684724 4.673478
## 426   -0.003923295 -3.063424 3.055578 -4.683026 4.675180
model5 <- arima(ts_inflation_change[311:730],order=c(2,0,0))
focst5 <- forecast(model5,h=6)
summary(focst5)
## 
## Forecast method: ARIMA(2,0,0) with non-zero mean
## 
## Model Information:
## 
## Call:
## arima(x = ts_inflation_change[311:730], order = c(2, 0, 0))
## 
## Coefficients:
##           ar1      ar2  intercept
##       -0.3196  -0.2865    -0.0050
## s.e.   0.0468   0.0467     0.0674
## 
## sigma^2 estimated as 4.906:  log likelihood = -930.07,  aic = 1868.15
## 
## Error measures:
##                      ME     RMSE      MAE      MPE     MAPE      MASE       ACF1
## Training set 0.00254465 2.214983 1.615356 53.48843 133.8038 0.5587509 -0.0643723
## 
## Forecasts:
##     Point Forecast     Lo 80    Hi 80     Lo 95    Hi 95
## 421     0.08070127 -2.757913 2.919316 -4.260585 4.421987
## 422    -0.44347184 -3.423576 2.536632 -5.001148 4.114204
## 423     0.11063824 -2.915064 3.136341 -4.516774 4.738051
## 424     0.08370913 -2.972010 3.139429 -4.589611 4.757029
## 425    -0.06645105 -3.122200 2.989298 -4.739816 4.606914
## 426    -0.01073694 -3.069111 3.047637 -4.688117 4.666643
model6 <- arima(ts_inflation_change[311:730],order=c(1,0,2))
focst6 <- forecast(model6,h=6)
summary(focst6)
## 
## Forecast method: ARIMA(1,0,2) with non-zero mean
## 
## Model Information:
## 
## Call:
## arima(x = ts_inflation_change[311:730], order = c(1, 0, 2))
## 
## Coefficients:
##          ar1      ma1      ma2  intercept
##       0.2398  -0.8080  -0.1920    -0.0047
## s.e.  0.1130   0.1141   0.1137     0.0012
## 
## sigma^2 estimated as 3.96:  log likelihood = -887.66,  aic = 1785.31
## 
## Error measures:
##                     ME     RMSE      MAE      MPE     MAPE      MASE         ACF1
## Training set -0.040883 1.990098 1.441103 44.32186 146.8271 0.4984767 -0.001788268
## 
## Forecasts:
##     Point Forecast     Lo 80    Hi 80     Lo 95    Hi 95
## 421   -0.650579993 -3.204009 1.902849 -4.555712 3.254552
## 422   -0.438510146 -3.372401 2.495381 -4.925510 4.048489
## 423   -0.108729360 -3.159273 2.941814 -4.774133 4.556674
## 424   -0.029655084 -3.086770 3.027460 -4.705109 4.645799
## 425   -0.010694788 -3.068187 3.046798 -4.686726 4.665336
## 426   -0.006148521 -3.063663 3.051365 -4.682213 4.669916
autoplot(ts(focst4[["mean"]]),main="6 month forecasts")+autolayer(ts(focst5[["mean"]]))+autolayer(ts(focst6[["mean"]]))

#########

model7 <- arima(ts_inflation_change[311:730],order=c(1,0,0))
focst7 <- forecast(model7,h=9)
summary(focst7)
## 
## Forecast method: ARIMA(1,0,0) with non-zero mean
## 
## Model Information:
## 
## Call:
## arima(x = ts_inflation_change[311:730], order = c(1, 0, 0))
## 
## Coefficients:
##           ar1  intercept
##       -0.2486    -0.0043
## s.e.   0.0473     0.0904
## 
## sigma^2 estimated as 5.347:  log likelihood = -948.07,  aic = 1902.13
## 
## Error measures:
##                       ME     RMSE     MAE      MPE     MAPE      MASE       ACF1
## Training set 0.001354483 2.312405 1.72315 81.28657 133.9653 0.5960368 -0.0710918
## 
## Forecasts:
##     Point Forecast     Lo 80    Hi 80     Lo 95    Hi 95
## 421   -0.360748069 -3.324214 2.602718 -4.892979 4.171483
## 422    0.084354337 -2.969300 3.138009 -4.585807 4.754516
## 423   -0.026289979 -3.085430 3.032850 -4.704841 4.652261
## 424    0.001214174 -3.058265 3.060693 -4.677855 4.680283
## 425   -0.005622855 -3.065123 3.053877 -4.684724 4.673478
## 426   -0.003923295 -3.063424 3.055578 -4.683026 4.675180
## 427   -0.004345775 -3.063847 3.055155 -4.683449 4.674758
## 428   -0.004240754 -3.063742 3.055260 -4.683344 4.674863
## 429   -0.004266860 -3.063768 3.055234 -4.683370 4.674836
model8 <- arima(ts_inflation_change[311:730],order=c(2,0,0))
focst8 <- forecast(model8,h=9)
summary(focst8)
## 
## Forecast method: ARIMA(2,0,0) with non-zero mean
## 
## Model Information:
## 
## Call:
## arima(x = ts_inflation_change[311:730], order = c(2, 0, 0))
## 
## Coefficients:
##           ar1      ar2  intercept
##       -0.3196  -0.2865    -0.0050
## s.e.   0.0468   0.0467     0.0674
## 
## sigma^2 estimated as 4.906:  log likelihood = -930.07,  aic = 1868.15
## 
## Error measures:
##                      ME     RMSE      MAE      MPE     MAPE      MASE       ACF1
## Training set 0.00254465 2.214983 1.615356 53.48843 133.8038 0.5587509 -0.0643723
## 
## Forecasts:
##     Point Forecast     Lo 80    Hi 80     Lo 95    Hi 95
## 421    0.080701275 -2.757913 2.919316 -4.260585 4.421987
## 422   -0.443471838 -3.423576 2.536632 -5.001148 4.114204
## 423    0.110638239 -2.915064 3.136341 -4.516774 4.738051
## 424    0.083709132 -2.972010 3.139429 -4.589611 4.757029
## 425   -0.066451046 -3.122200 2.989298 -4.739816 4.606914
## 426   -0.010736943 -3.069111 3.047637 -4.688117 4.666643
## 427    0.014479316 -3.044115 3.073073 -4.663237 4.692195
## 428   -0.009544613 -3.068237 3.049148 -4.687412 4.668322
## 429   -0.009090611 -3.067839 3.049657 -4.687042 4.668861
model9 <- arima(ts_inflation_change[311:730],order=c(1,0,2))
focst9 <- forecast(model9,h=9)
summary(focst9)
## 
## Forecast method: ARIMA(1,0,2) with non-zero mean
## 
## Model Information:
## 
## Call:
## arima(x = ts_inflation_change[311:730], order = c(1, 0, 2))
## 
## Coefficients:
##          ar1      ma1      ma2  intercept
##       0.2398  -0.8080  -0.1920    -0.0047
## s.e.  0.1130   0.1141   0.1137     0.0012
## 
## sigma^2 estimated as 3.96:  log likelihood = -887.66,  aic = 1785.31
## 
## Error measures:
##                     ME     RMSE      MAE      MPE     MAPE      MASE         ACF1
## Training set -0.040883 1.990098 1.441103 44.32186 146.8271 0.4984767 -0.001788268
## 
## Forecasts:
##     Point Forecast     Lo 80    Hi 80     Lo 95    Hi 95
## 421   -0.650579993 -3.204009 1.902849 -4.555712 3.254552
## 422   -0.438510146 -3.372401 2.495381 -4.925510 4.048489
## 423   -0.108729360 -3.159273 2.941814 -4.774133 4.556674
## 424   -0.029655084 -3.086770 3.027460 -4.705109 4.645799
## 425   -0.010694788 -3.068187 3.046798 -4.686726 4.665336
## 426   -0.006148521 -3.063663 3.051365 -4.682213 4.669916
## 427   -0.005058424 -3.062574 3.052457 -4.681124 4.671008
## 428   -0.004797043 -3.062312 3.052718 -4.680863 4.671269
## 429   -0.004734369 -3.062250 3.052781 -4.680801 4.671332
autoplot(ts(focst7[["mean"]]),main="9 month forecasts")+autolayer(ts(focst8[["mean"]]))+autolayer(ts(focst9[["mean"]]))


As the forecast horizon increases, all three of the models forecasts seem to converge, and so it is important that we select the optimal forecast horizon and the optimal model in order to get meaningful results.