setwd("C:/Users/victo/Downloads/ECON 4310/Paper")We are using various R packages that need to be installed and added. You might have to uncomment/run the chunk below. Alternatively, you may be able to install the missing packages through RStudio, selecting “packages” and “install”, entering the names.
#install.packages("invgamma")
library(invgamma)
#install.packages("matlib")
library(matlib)
#install.packages("MASS")
library(MASS)
#install.packages("MCMCpack")
library(MCMCpack) #for Inverse Wishart## Warning: package 'MCMCpack' was built under R version 4.2.3## Loading required package: coda## Warning: package 'coda' was built under R version 4.2.3## ##
## ## Markov Chain Monte Carlo Package (MCMCpack)## ## Copyright (C) 2003-2023 Andrew D. Martin, Kevin M. Quinn, and Jong Hee Park## ##
## ## Support provided by the U.S. National Science Foundation## ## (Grants SES-0350646 and SES-0350613)
## ####
## Attaching package: 'MCMCpack'## The following objects are masked from 'package:invgamma':
##
## dinvgamma, rinvgamma#install.packages("mvtnorm")
library(mvtnorm) #for multivariate Normal
#install.packages("HypergeoMat")
library(HypergeoMat) #for multivariate gamma function (VAR-MDD calculation)## Warning: package 'HypergeoMat' was built under R version 4.2.3#install.packages("vars")
library(vars)## Warning: package 'vars' was built under R version 4.2.3## Loading required package: strucchange## Warning: package 'strucchange' was built under R version 4.2.3## Loading required package: zoo## Warning: package 'zoo' was built under R version 4.2.3##
## Attaching package: 'zoo'## The following objects are masked from 'package:base':
##
## as.Date, as.Date.numeric## Loading required package: sandwich## Warning: package 'sandwich' was built under R version 4.2.3## Loading required package: urca## Warning: package 'urca' was built under R version 4.2.3## Loading required package: lmtest## Warning: package 'lmtest' was built under R version 4.2.3Load monthly data and trim to 1997, Q2 (max MPR data):
data <- read.csv("Chile_1986_Monthly.csv", header=TRUE)
# head(data)
#na.omit(as.numeric(data$CPI))
#na.omit(data$CPI)
UE <- ts(data$UE, frequency=12,start=c(1986,1))
CPI <- ts(data$CPI, frequency=12, start=c(1986,1))
UE <- window(UE,start=c(1997,2),end=c(2022,12))
CPI <- window(CPI,start=c(1997,2),end=c(2022,12))Convert monthly to quarterly data
MonthlyToQuarterly <- function(MonthlyData,start,namevec){
## Assume that first obs is first month of quarter
## Assume that last obs is last month of quarter
nobsMonth <- length(MonthlyData)
nobsQuart <- nobsMonth/3
firstyear <- start[1]
firstquarter <- ((start[2]-1) %/% 3)+1
QuarterlyData <- matrix(0,nrow=nobsQuart,ncol=1)
for(i in 1:nobsQuart){
QuarterlyData[i] <- mean(MonthlyData[(1+3*(i-1)):(3*i)],na.rm=FALSE)
}
QuarterlyData <- ts(QuarterlyData,frequency=4,start=c(firstyear,firstquarter),
names=namevec)
return(QuarterlyData)
}
UEq <- MonthlyToQuarterly(UE,c(1996,1),"UE")
UEq <- window(UEq,start=c(1997,2))
#print(UEq)
CPIq <- MonthlyToQuarterly(CPI,c(1996,1),"CPI")
CPIq <- window(CPIq,start=c(1997,2))
#print(CPIq)Load Monetary Policy Rate data (MPR)
data <- read.csv("Chile_MPR_monthly_F1997.csv", header = T)
MPR <- ts(data$MPR, frequency = 12, start = c(1997, 2))
MPRq <- MonthlyToQuarterly(MPR, c(1997,2), "MPR")
MPRq <- window(MPRq, start=c(1997,2), end=c(2021,3))
#print(MPRq)Load the GDP data (unnecessary for monetary policy shock)
data <- read.csv("Chile_1996_Quarterly.csv", header=TRUE)
GDP <- gsub(",","",data$GDP)
GDP <- as.numeric(GDP)
GDP <- ts(GDP, frequency=4, start=c(1997,2))
#head(GDP)
GDPlog <- log(GDP)
GDPlogdiff <- diff(GDPlog)
GDPgrowth <- 400*GDPlogdiff
GDPgrowth <- ts(GDPgrowth, frequency = 4,start=c(1997,2))# Prepare data
Tall <- length(UEq) # total number of observations, starting 1996:Q1
T1 <- (1+4) # estimation sample will start in T1 = 1997:Q1
p <- 2 # number of lags
n <- 3
yyall <- cbind(UEq,CPIq,MPRq)
yy <- yyall[T1:Tall,]
xx <- matrix(1,nrow=(Tall-T1+1), ncol=(n*p+1))
for(i in 1:p){
# first column of xx is vector of 1s
# followed by y(t-1), y(t-2), ..., y(t-p)
xx[,(1+(i-1)*n+1):(1+i*n)] <- yyall[(T1-i):(Tall-i),]
}
head(yy)## UEq CPIq MPRq
## [1,] 10.779787 0.1090467 8.500000
## [2,] 10.851437 0.2595300 10.706667
## [3,] 9.640314 0.2023733 8.610000
## [4,] 9.161562 0.5886367 6.966667
## [5,] 9.878699 0.1912600 5.550000
## [6,] 10.289450 0.4904367 5.000000head(xx)## [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,] 1 9.386524 0.3611800 8.450000 8.214573 0.09094667 6.633333
## [2,] 1 10.779787 0.1090467 8.500000 9.386524 0.36118000 8.450000
## [3,] 1 10.851437 0.2595300 10.706667 10.779787 0.10904667 8.500000
## [4,] 1 9.640314 0.2023733 8.610000 10.851437 0.25953000 10.706667
## [5,] 1 9.161562 0.5886367 6.966667 9.640314 0.20237333 8.610000
## [6,] 1 9.878699 0.1912600 5.550000 9.161562 0.58863667 6.966667Compute OLS/ML estimators of Phi and Sigma
Phi_hat <- solve(t(xx)%*%xx,t(xx)%*%yy)
Sigma_hat <- (t(yy)%*%yy - t(yy)%*%xx%*%Phi_hat)/(Tall-T1+1)
print(Phi_hat)## UEq CPIq MPRq
## [1,] 0.713705868 0.34727630 0.985396748
## [2,] 1.281559771 -0.03755187 0.160662281
## [3,] -0.150419435 0.21643396 0.043414486
## [4,] 0.016628894 -0.06506056 1.257598589
## [5,] -0.385509465 0.02874264 -0.220691200
## [6,] 0.343584324 0.13338757 0.004729429
## [7,] 0.003023265 0.04488955 -0.399775258print(Sigma_hat)## UEq CPIq MPRq
## UEq 0.25850085 -0.0138286788 0.0300430397
## CPIq -0.01382868 0.0956714142 0.0007385997
## MPRq 0.03004304 0.0007385997 0.4972743691Cholesky Factorization in R
Sigma_hat_tr <- t(chol(Sigma_hat))
print(Sigma_hat_tr)## UEq CPIq MPRq
## UEq 0.50842979 0.000000000 0.0000000
## CPIq -0.02719880 0.308109785 0.0000000
## MPRq 0.05908985 0.007613431 0.7026555# Verify that Sigma_hat_tr*t(Sigma_hat_tr) = Sigma_hat
#print(Sigma_hat)
#print(Sigma_hat_tr %*% t(Sigma_hat_tr))VAR Model in level
nhor = 12
chile <- cbind(UEq, CPIq, MPRq)
chile.lev <- chile
# lag length
VARselect(chile.lev, lag.max = 4, type = "const", season = 4)## $selection
## AIC(n) HQ(n) SC(n) FPE(n)
## 2 2 2 2
##
## $criteria
## 1 2 3 4
## AIC(n) -3.92234275 -4.1933284 -4.06962271 -4.10633194
## HQ(n) -3.69283845 -3.8654651 -3.64340044 -3.58175070
## SC(n) -3.35415987 -3.3816386 -3.01442593 -2.80762823
## FPE(n) 0.01981108 0.0151326 0.01717551 0.01663421#estimation
var.model_lev <- VAR(chile.lev, p = 2, type = "const", season = 4)
# forecast of lev data
var.pred <- predict(var.model_lev, n.ahead = nhor)
x11(); par(mai=rep(0.4, 4)); plot(var.pred)
x11(); par(mai=rep(0.4, 4)); fanchart(var.pred)