Final Paper: VAR by OLS (IRF)
Set work directory
setwd("C:/Users/victo/Downloads/ECON 4310/Paper")Loading Packages
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.3Load Data
Load 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))Estimate VAR
# Prepare data
Tall <- length(UEq) # total number of observations, starting 1997:Q2
T1 <- (1+3) # estimation sample will start in T1 = 1998:Q1
p <- 2 # number of lags
n <- 3 # number of series in VAR
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,] 9.386524 0.3611800 8.450000
## [2,] 10.779787 0.1090467 8.500000
## [3,] 10.851437 0.2595300 10.706667
## [4,] 9.640314 0.2023733 8.610000
## [5,] 9.161562 0.5886367 6.966667
## [6,] 9.878699 0.1912600 5.550000head(xx)## [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,] 1 8.214573 0.09094667 6.633333 6.791117 0.55271333 6.593333
## [2,] 1 9.386524 0.36118000 8.450000 8.214573 0.09094667 6.633333
## [3,] 1 10.779787 0.10904667 8.500000 9.386524 0.36118000 8.450000
## [4,] 1 10.851437 0.25953000 10.706667 10.779787 0.10904667 8.500000
## [5,] 1 9.640314 0.20237333 8.610000 10.851437 0.25953000 10.706667
## [6,] 1 9.161562 0.58863667 6.966667 9.640314 0.20237333 8.610000Compute 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.720572571 0.34982157 1.01441247
## [2,] 1.298268331 -0.03135854 0.23126546
## [3,] -0.157578608 0.21378029 0.01316291
## [4,] 0.015456130 -0.06549527 1.25264299
## [5,] -0.404090542 0.02185523 -0.29920682
## [6,] 0.361527326 0.14003846 0.08054883
## [7,] 0.006729597 0.04626337 -0.38411389print(Sigma_hat)## UEq CPIq MPRq
## UEq 0.25776658 -0.012946675 0.038122113
## CPIq -0.01294668 0.094937320 0.003842696
## MPRq 0.03812211 0.003842696 0.527514891Cholesky Factorization in R
Sigma_hat_tr <- t(chol(Sigma_hat))
print(Sigma_hat_tr)## UEq CPIq MPRq
## UEq 0.50770718 0.00000000 0.0000000
## CPIq -0.02550028 0.30706197 0.0000000
## MPRq 0.07508681 0.01875006 0.7221671# Verify that Sigma_hat_tr*t(Sigma_hat_tr) = Sigma_hat
#print(Sigma_hat)
#print(Sigma_hat_tr %*% t(Sigma_hat_tr))Procedure to Compute IRF to Shock sh_ind
IRFVARp <- function(Phi,Sigma_tr,sh_ind,Hmax){
n <- ncol(Phi)
k <- nrow(Phi)
p <- k/n
yy <- matrix(0, nrow=(Hmax+1), ncol=n )
# define xx without intercept
xt <- matrix(0, nrow=1, ncol=n*p)
# Impact effect
yy[1,] <- t(Sigma_tr[,sh_ind])
# loop to generate the impulse responses
xt[,1:n] <- yy[1,]
for(t in 2:(Hmax+1)){
yy[t,] <- xt%*%Phi
# update the x(t) vector -> x(t+1)
if(p > 1){
xt[,(n+1):(p*n)] <- xt[,1:((p-1)*n)]
xt[,1:n] <- yy[t,]}
else{
xt[,1:n] <- yy[t,]}
}
return(yy)
}
Hmax <- 120
sh_ind <- 1
# we are dropping the intercept when computing the IRF
VARirf_sh1 <- IRFVARp(Phi_hat[2:nrow(Phi_hat),],Sigma_hat_tr,sh_ind,Hmax)
sh_ind <- 2
# we are dropping the intercept when computing the IRF
VARirf_sh2 <- IRFVARp(Phi_hat[2:nrow(Phi_hat),],Sigma_hat_tr,sh_ind,Hmax)
plot(0:Hmax,VARirf_sh1[,1],type="l",cex.axis=1.5,cex.lab=1.5,xlab="Horizon h", ylab="Unemployment Response",col="blue",lwd=2)
abline(h=0,col="black",lwd=1)
plot(0:Hmax,VARirf_sh1[,2],type="l",cex.axis=1.5,cex.lab=1.5,xlab="Horizon h", ylab="CPI Response", col="blue", lwd=2)
abline(h=0,col="black",lwd=1)
plot(0:Hmax,VARirf_sh1[,3],type="l",cex.axis=1.5,cex.lab=1.5,xlab="Horizon h", ylab="MPR Response", col="blue", lwd=2)
abline(h=0,col="black",lwd=1)