library(forecast)
## Registered S3 method overwritten by 'quantmod':
## method from
## as.zoo.data.frame zoo
library(lmtest)
## Loading required package: zoo
##
## Attaching package: 'zoo'
## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric
library(car)
## Loading required package: carData
library(vars)
## Loading required package: MASS
## Loading required package: strucchange
## Loading required package: sandwich
## Loading required package: urca
library(dygraphs)
order_checker <- function(Y, max.order = 3, xreg = NULL, I = 0, include.drift = FALSE){
params <- list()
for (j in c("ar", "ma")){
for (i in 1:max.order){
params[[paste0(j, i)]] <- c(0, NA)
}
}
if (!is.null(xreg)){
for (j in colnames(xreg)){
params[[j]] <- NA
}
}
params$intercept <- NA
combinations <- expand.grid(params)
res <- data.frame(BIC = rep(NA, nrow(combinations)), AIC = NA)
pb <- txtProgressBar(1, nrow(combinations), style = 3)
for (i in 1:nrow(combinations)){
curr_mod <- try(Arima(Y, order = c(max.order, I, max.order), fixed = c(combinations[i, ]),
xreg = xreg, method = "ML", include.drift = include.drift), silent = TRUE)
if ("try-error" %in% attr(curr_mod, "class")){
res$BIC[i] <- NA
res$AIC[i] <- NA
} else {
res$BIC[i] <- BIC(curr_mod)
res$AIC[i] <- AIC(curr_mod)
}
setTxtProgressBar(pb, i)
}
res <- merge(res, combinations, by = 0)
res$Row.names <- NULL
return(res)
}
inflation=read.csv("inflation.csv")
money=read.csv("Money.csv")
inflation_ts=ts(inflation[,c("s1")],start=c(1997,01),frequency=12)#HCIP annual change
money_ts=ts(money[,c("s1")],start=c(1981,01),end=c(2021,03),frequency=12)
money_ts2=window(money_ts,start=c(1997,01),end=c(2021,03))
inflation_ts2=window(inflation_ts,start=c(1997,01),end=c(2019,12))
res <- order_checker(inflation_ts2, 4)
##
|
| | 0%
|
|= | 1%
|
|= | 2%
|
|== | 2%
|
|== | 3%
|
|== | 4%
|
|=== | 4%
|
|=== | 5%
|
|==== | 5%
|
|==== | 6%
|
|===== | 7%
|
|===== | 8%
|
|====== | 8%
|
|====== | 9%
|
|======= | 9%
|
|======= | 10%
|
|======= | 11%
|
|======== | 11%
|
|======== | 12%
|
|========= | 12%
|
|========= | 13%
|
|========== | 14%
|
|========== | 15%
|
|=========== | 15%
|
|=========== | 16%
|
|============ | 16%
|
|============ | 17%
|
|============ | 18%
|
|============= | 18%
|
|============= | 19%
|
|============== | 20%
|
|=============== | 21%
|
|=============== | 22%
|
|================ | 22%
|
|================ | 23%
|
|================ | 24%
|
|================= | 24%
|
|================= | 25%
|
|================== | 25%
|
|================== | 26%
|
|=================== | 27%
|
|=================== | 28%
|
|==================== | 28%
|
|==================== | 29%
|
|===================== | 29%
|
|===================== | 30%
|
|===================== | 31%
|
|====================== | 31%
|
|====================== | 32%
|
|======================= | 32%
|
|======================= | 33%
|
|======================== | 34%
|
|======================== | 35%
|
|========================= | 35%
|
|========================= | 36%
|
|========================== | 36%
|
|========================== | 37%
|
|========================== | 38%
|
|=========================== | 38%
|
|=========================== | 39%
|
|============================ | 40%
|
|============================= | 41%
|
|============================= | 42%
|
|============================== | 42%
|
|============================== | 43%
|
|============================== | 44%
|
|=============================== | 44%
|
|=============================== | 45%
|
|================================ | 45%
|
|================================ | 46%
|
|================================= | 47%
|
|================================= | 48%
|
|================================== | 48%
|
|================================== | 49%
|
|=================================== | 49%
|
|=================================== | 50%
|
|=================================== | 51%
|
|==================================== | 51%
|
|==================================== | 52%
|
|===================================== | 52%
|
|===================================== | 53%
|
|====================================== | 54%
|
|====================================== | 55%
|
|======================================= | 55%
|
|======================================= | 56%
|
|======================================== | 56%
|
|======================================== | 57%
|
|======================================== | 58%
|
|========================================= | 58%
|
|========================================= | 59%
|
|========================================== | 60%
|
|=========================================== | 61%
|
|=========================================== | 62%
|
|============================================ | 62%
|
|============================================ | 63%
|
|============================================ | 64%
|
|============================================= | 64%
|
|============================================= | 65%
|
|============================================== | 65%
|
|============================================== | 66%
|
|=============================================== | 67%
|
|=============================================== | 68%
|
|================================================ | 68%
|
|================================================ | 69%
|
|================================================= | 69%
|
|================================================= | 70%
|
|================================================= | 71%
|
|================================================== | 71%
|
|================================================== | 72%
|
|=================================================== | 72%
|
|=================================================== | 73%
|
|==================================================== | 74%
|
|==================================================== | 75%
|
|===================================================== | 75%
|
|===================================================== | 76%
|
|====================================================== | 76%
|
|====================================================== | 77%
|
|====================================================== | 78%
|
|======================================================= | 78%
|
|======================================================= | 79%
|
|======================================================== | 80%
|
|========================================================= | 81%
|
|========================================================= | 82%
|
|========================================================== | 82%
|
|========================================================== | 83%
|
|========================================================== | 84%
|
|=========================================================== | 84%
|
|=========================================================== | 85%
|
|============================================================ | 85%
|
|============================================================ | 86%
|
|============================================================= | 87%
|
|============================================================= | 88%
|
|============================================================== | 88%
|
|============================================================== | 89%
|
|=============================================================== | 89%
|
|=============================================================== | 90%
|
|=============================================================== | 91%
|
|================================================================ | 91%
|
|================================================================ | 92%
|
|================================================================= | 92%
|
|================================================================= | 93%
|
|================================================================== | 94%
|
|================================================================== | 95%
|
|=================================================================== | 95%
|
|=================================================================== | 96%
|
|==================================================================== | 96%
|
|==================================================================== | 97%
|
|==================================================================== | 98%
|
|===================================================================== | 98%
|
|===================================================================== | 99%
|
|======================================================================| 100%
final_us_mod <- Arima(inflation_ts2, order = c(4, 0, 4), fixed = res[which.min(res$BIC), 3:ncol(res)])
Acf(inflation_ts2)
mod=Arima(inflation_ts2,c(1,0,1),seasonal=c(0,0,0))
Acf(residuals(mod),main="ARIMA(1,0,1)")
mod_012s=Arima(inflation_ts2,c(0,1,2),seasonal=c(0,0,1))
mod_102s=Arima(inflation_ts2,c(1,0,2),seasonal=c(0,0,1))
mod_101s=Arima(inflation_ts2,c(1,0,1),seasonal=list(order=c(0,0,1),period=12))
mod_1011s=Arima(inflation_ts2,c(1,0,1),seasonal=list(order=c(1,0,1),period=12))
mod_101s4=Arima(inflation_ts2,c(1,0,1),seasonal=list(order=c(0,0,1),period=4))
mod_404s=Arima(inflation_ts2,c(4,0,4),seasonal=list(order=c(0,0,1),period=12),fixed =c(0,NA,NA,NA,NA,NA,NA,0,NA,0))
summary(mod_101s)
## Series: inflation_ts2
## ARIMA(1,0,1)(0,0,1)[12] with non-zero mean
##
## Coefficients:
## ar1 ma1 sma1 mean
## 0.9780 0.1259 -0.5320 1.6323
## s.e. 0.0119 0.0552 0.0482 0.2789
##
## sigma^2 estimated as 0.04599: log likelihood=32.17
## AIC=-54.34 AICc=-54.12 BIC=-36.24
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -0.0009147759 0.2128926 0.1633434 NaN Inf 0.2047609 0.01675587
mod1=Arima(inflation_ts2,c(1,0,1))
plot(residuals(mod1))
criteriatab=cbind(AIC(mod1,mod_102s,mod_101s,mod_101s4,final_us_mod,mod_404s,mod_1011s),
BIC(mod1,mod_102s,mod_101s,mod_101s4,final_us_mod,mod_404s,mod_1011s)[,2])%>%
dplyr::rename(BIC=starts_with("BIC"))
criteriatab
## df AIC BIC
## mod1 4 15.83611 30.31771
## mod_102s 6 -55.06536 -33.34296
## mod_101s 5 -54.34189 -36.23989
## mod_101s4 5 15.47578 33.57778
## final_us_mod 8 -18.35517 10.60804
## mod_404s 8 -57.90739 -28.94418
## mod_1011s 6 -53.23025 -31.50784
par(mfrow=c(1,1))
Acf(residuals(mod_101s),main="SARIMA(1,0,1)(0,0,1)12",ylim=c(-0.4,0.2))
Pacf(residuals(mod_101s),main="SARIMA(1,0,1)(0,0,1)12",ylim=c(-0.4,0.2))
plot(residuals(mod_101s))
table <- data.frame(c(BIC(mod),BIC(mod_101s)),c(AIC(mod), AIC(mod_101s)))
colnames(table) <- c("BIC", "AIC")
rownames(table) <- c("ARIMA(1,0,1)","SARIMA(1,0,1)(0,0,1)[12]")
table
## BIC AIC
## ARIMA(1,0,1) 30.31771 15.83611
## SARIMA(1,0,1)(0,0,1)[12] -36.23989 -54.34189
plot(inflation_ts2)
lines(fitted(mod_101s),col="red")
plot(forecast(mod_101s,20))
#lines(test,col="red")
test=window(inflation_ts,start=c(2019,12))
##################################################### with covid ###############################
inflation_ts2=window(inflation_ts,start=c(1997,01),end=c(2019,12))
test=window(inflation_ts,start=c(2019,12),end=c(2021,4))
mod_101s=Arima(inflation_ts2,c(1,0,1),seasonal=list(order=c(0,0,1),period=12))
summary(mod_101s)
## Series: inflation_ts2
## ARIMA(1,0,1)(0,0,1)[12] with non-zero mean
##
## Coefficients:
## ar1 ma1 sma1 mean
## 0.9780 0.1259 -0.5320 1.6323
## s.e. 0.0119 0.0552 0.0482 0.2789
##
## sigma^2 estimated as 0.04599: log likelihood=32.17
## AIC=-54.34 AICc=-54.12 BIC=-36.24
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -0.0009147759 0.2128926 0.1633434 NaN Inf 0.2047609 0.01675587
Acf(residuals(mod_101s))
Pacf(residuals(mod_101s))
plot(residuals(mod_101s))
plot(inflation_ts2)
lines(fitted(mod_101s),col="red")
#plot(forecast(mod_101s,20))
#lines(test,col="red")
forc=forecast(mod_101s,16)
dygraph(inflation_ts)
#dygraph(forecast(mod_101s,20))
gen_array <- function(forecast_obj){
actuals <- forecast_obj$x
lower <- forecast_obj$lower[,2]
upper <- forecast_obj$upper[,2]
point_forecast <- forecast_obj$mean
cbind(actuals, lower, upper, point_forecast)
}
forc=gen_array(forc)
dygraph(test)
df=as.data.frame(test)
ntest=c(rep(NA,times=nrow(forc)-nrow(df)),as.vector(df$x))
ntest=ts(ntest,start=c(1997,01),end=c(2021,4),frequency=12)
cast=as.data.frame(forc)
cast$test=ntest
forecast=ts(cast, start=c(1997,01),frequency=12)
forecast[[276,4]]=1.3
dygraph(forecast, main = "Forecast 2020",ylab="Percentage Change") %>%
dyRangeSelector() %>%
dyRangeSelector(height = 40,
dateWindow = c("2011-04-01", "2019-4-01")) %>%
dySeries(name = "actuals", label = "Actual") %>%
#dySeries(name="ntest", label = "Predicted") %>%
dySeries(c("lower","point_forecast","upper"), label = "Predicted") %>%
dySeries(name="test", label = "Test Data") %>%
dyLegend(show = "always", hideOnMouseOut = FALSE) %>%
dyHighlight(highlightCircleSize = 5,
highlightSeriesOpts = list(strokeWidth = 2)) %>%
dyOptions(axisLineColor = "navy", gridLineColor = "grey")
#################################################### 2019 Onward ####################################
inflation_ts2=window(inflation_ts,start=c(1997,01),end=c(2018,01))
test=window(inflation_ts,start=c(2017,12),end=c(2019,01))
mod_101s=Arima(inflation_ts2,c(1,0,1),seasonal=list(order=c(0,0,1),period=12))
summary(mod_101s)
## Series: inflation_ts2
## ARIMA(1,0,1)(0,0,1)[12] with non-zero mean
##
## Coefficients:
## ar1 ma1 sma1 mean
## 0.9778 0.1468 -0.5135 1.6398
## s.e. 0.0125 0.0580 0.0532 0.3040
##
## sigma^2 estimated as 0.04618: log likelihood=29
## AIC=-47.99 AICc=-47.75 BIC=-30.33
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE ACF1
## Training set -0.0009729089 0.2131853 0.1639349 NaN Inf 0.2022955 0.01821875
Acf(residuals(mod_101s))
Pacf(residuals(mod_101s))
plot(residuals(mod_101s))
plot(inflation_ts2)
lines(fitted(mod_101s),col="red")
forc=forecast(mod_101s,13)
gen_array <- function(forecast_obj){
actuals <- forecast_obj$x
lower <- forecast_obj$lower[,2]
upper <- forecast_obj$upper[,2]
point_forecast <- forecast_obj$mean
cbind(actuals, lower, upper, point_forecast)
}
forc=gen_array(forc)
df=as.data.frame(test)
ntest=c(rep(NA,times=nrow(forc)-nrow(df)),as.vector(df$x))
ntest=ts(ntest,start=c(1997,01),end=c(2019,02),frequency=12)
cast=as.data.frame(forc)
cast$test=ntest
forecast=ts(cast, start=c(1997,01),frequency=12)
forecast[(nrow(forecast)-13),4]=1.3
dygraph(forecast, main = "Forecast 2019",ylab="Percentage Change") %>%
dyRangeSelector() %>%
dyRangeSelector(height = 40,
dateWindow = c("2016-01-01", "2019-01-01")) %>%
dySeries(name = "actuals", label = "Actual") %>%
#dySeries(name="test", label = "Predicted") %>%
dySeries(c("lower","point_forecast","upper"), label = "Predicted") %>%
dySeries(name="test", label = "Test Data") %>%
dyLegend(show = "always", hideOnMouseOut = FALSE) %>%
dyHighlight(highlightCircleSize = 5,
highlightSeriesOpts = list(strokeWidth = 2)) %>%
dyOptions(axisLineColor = "navy", gridLineColor = "grey")
######################################################### 2005#####################
inflation_ts2=window(inflation_ts,start=c(1997,01),end=c(2004,12))
test=window(inflation_ts,start=c(2004,12),end=c(2006,01))
mod_101s=Arima(inflation_ts2,c(1,0,1),seasonal=list(order=c(0,0,1),period=12))
summary(mod_101s)
## Series: inflation_ts2
## ARIMA(1,0,1)(0,0,1)[12] with non-zero mean
##
## Coefficients:
## ar1 ma1 sma1 mean
## 0.9557 0.1999 -0.5242 1.8782
## s.e. 0.0306 0.1126 0.1023 0.2096
##
## sigma^2 estimated as 0.03107: log likelihood=29.48
## AIC=-48.95 AICc=-48.28 BIC=-36.13
##
## Training set error measures:
## ME RMSE MAE MPE MAPE MASE
## Training set -0.0002251893 0.1725538 0.1419162 -1.393153 8.319214 0.2987709
## ACF1
## Training set 0.00121644
Acf(residuals(mod_101s))
Pacf(residuals(mod_101s))
plot(residuals(mod_101s))
plot(inflation_ts2)
lines(fitted(mod_101s),col="red")
forc=forecast(mod_101s,13)
plot(forecast(mod_101s,10))
lines(test,col="red")
gen_array <- function(forecast_obj){
actuals <- forecast_obj$x
lower <- forecast_obj$lower[,2]
upper <- forecast_obj$upper[,2]
point_forecast <- forecast_obj$mean
cbind(actuals, lower, upper, point_forecast)
}
forc=gen_array(forc)
dygraph(test)
df=as.data.frame(test)
ntest=c(rep(NA,times=nrow(forc)-nrow(df)),as.vector(df$x))
ntest=ts(ntest,start=c(1997,01),end=c(2006,01),frequency=12)
cast=as.data.frame(forc)
cast$test=ntest
forecast=ts(cast, start=c(1997,01),frequency=12)
forecast[(nrow(forecast)-13),4]=2.3
dygraph(forecast, main = "Forecast 2005",ylab="Percentage Change") %>%
dyRangeSelector() %>%
dyRangeSelector(height = 40,
dateWindow = c("2003-04-01", "2005-01-01")) %>%
dySeries(name = "actuals", label = "Actual") %>%
#dySeries(name="test", label = "Predicted") %>%
dySeries(c("lower","point_forecast","upper"), label = "Predicted") %>%
dySeries(name="test", label = "Test Data") %>%
dyLegend(show = "always", hideOnMouseOut = FALSE) %>%
dyHighlight(highlightCircleSize = 5,
highlightSeriesOpts = list(strokeWidth = 2)) %>%
dyOptions(axisLineColor = "navy", gridLineColor = "grey")
#######################################mutlivariate############################################################
data_ts = window(ts.union( money_ts,inflation_ts),start=c(1997,01),end=c(2021,03))
var_1=VAR(data_ts,p=1)
summary(var_1)
##
## VAR Estimation Results:
## =========================
## Endogenous variables: money_ts, inflation_ts
## Deterministic variables: const
## Sample size: 290
## Log Likelihood: -213.935
## Roots of the characteristic polynomial:
## 0.989 0.9579
## Call:
## VAR(y = data_ts, p = 1)
##
##
## Estimation results for equation money_ts:
## =========================================
## money_ts = money_ts.l1 + inflation_ts.l1 + const
##
## Estimate Std. Error t value Pr(>|t|)
## money_ts.l1 0.99063 0.01015 97.567 <2e-16 ***
## inflation_ts.l1 -0.01001 0.03153 -0.318 0.751
## const 0.08800 0.07220 1.219 0.224
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Residual standard error: 0.4811 on 287 degrees of freedom
## Multiple R-Squared: 0.9718, Adjusted R-squared: 0.9716
## F-statistic: 4952 on 2 and 287 DF, p-value: < 2.2e-16
##
##
## Estimation results for equation inflation_ts:
## =============================================
## inflation_ts = money_ts.l1 + inflation_ts.l1 + const
##
## Estimate Std. Error t value Pr(>|t|)
## money_ts.l1 0.005273 0.005443 0.969 0.333
## inflation_ts.l1 0.956287 0.016901 56.582 <2e-16 ***
## const 0.039325 0.038703 1.016 0.310
## ---
## Signif. codes: 0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
##
##
## Residual standard error: 0.2579 on 287 degrees of freedom
## Multiple R-Squared: 0.9213, Adjusted R-squared: 0.9207
## F-statistic: 1679 on 2 and 287 DF, p-value: < 2.2e-16
##
##
##
## Covariance matrix of residuals:
## money_ts inflation_ts
## money_ts 0.231482 -0.009544
## inflation_ts -0.009544 0.066513
##
## Correlation matrix of residuals:
## money_ts inflation_ts
## money_ts 1.00000 -0.07692
## inflation_ts -0.07692 1.00000
plot(residuals(var_1)[,1],type="l")
plot(residuals(var_1)[,2],type="l")
################### plots ##############################
m=dygraph(money_ts, main = "Monetary Supply", ylab = "Percentage Change") %>%
dyRangeSelector() %>%
dyRangeSelector(height = 40,
dateWindow = c("2011-04-01", "2021-03-31")) %>%
dySeries(c("V1"), label = "Supply") %>%
dyLegend(show = "always", hideOnMouseOut = FALSE) %>%
dyHighlight(highlightCircleSize = 5,
highlightSeriesOpts = list(strokeWidth = 2)) %>%
dyOptions(axisLineColor = "navy", gridLineColor = "grey")
i=dygraph(inflation_ts, main = "Inflation", ylab = "Percentage Change") %>%
dyRangeSelector() %>%
dyRangeSelector(height = 40,
dateWindow = c("2011-04-01", "2021-03-31")) %>%
dySeries(c("V1"), label = "Inflation") %>%
dyLegend(show = "always", hideOnMouseOut = FALSE) %>%
dyHighlight(highlightCircleSize = 5,
highlightSeriesOpts = list(strokeWidth = 2)) %>%
dyOptions(axisLineColor = "navy", gridLineColor = "grey")
dy_graph <- list(
i,
m
) # end list
# render the dygraphs objects using htmltools
htmltools::browsable(htmltools::tagList(dy_graph))