Final Paper: VAR by OLS
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 1996 (max quarterly 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(1996,1),end=c(2022,12))
CPI <- window(CPI,start=c(1996,1),end=c(2022,12))
#as.numeric(CPI)
#na.omit(CPI)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")
#UE <- window(UEq,start=c(1996,1),end=c(2022,4))
CPIq <- MonthlyToQuarterly(CPI,c(1996,1),"CPI")
#CPI <- window(CPIq,start=c(1996,1), end=c(2022,4))Load the quarterly data
data <- read.csv("Chile_1996_Quarterly.csv", header=TRUE)
GDP <- gsub(",","",data$GDP)
GDP <- as.numeric(GDP)
GDP <- ts(GDP, frequency=4, start=c(1996,1))
#head(GDP)
GDPlog <- log(GDP)
GDPlogdiff <- diff(GDPlog)
GDPgrowth <- 400*GDPlogdiff
GDPgrowth <- ts(GDPgrowth, frequency = 4,start=c(1996,1))Estimate VAR
# 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,GDP)
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 GDP
## [1,] 6.070701 0.5512667 8959.62
## [2,] 6.290691 0.2665900 8830.08
## [3,] 6.182505 0.6460100 8534.64
## [4,] 5.864666 0.4996967 9622.65
## [5,] 5.521886 0.3208600 9659.35
## [6,] 5.815567 0.3059367 9497.54head(xx)## [,1] [,2] [,3] [,4] [,5] [,6] [,7]
## [1,] 1 5.921170 0.4910567 8624.09 6.236843 0.4023400 7640.74
## [2,] 1 6.070701 0.5512667 8959.62 5.921170 0.4910567 8624.09
## [3,] 1 6.290691 0.2665900 8830.08 6.070701 0.5512667 8959.62
## [4,] 1 6.182505 0.6460100 8534.64 6.290691 0.2665900 8830.08
## [5,] 1 5.864666 0.4996967 9622.65 6.182505 0.6460100 8534.64
## [6,] 1 5.521886 0.3208600 9659.35 5.864666 0.4996967 9622.65Compute 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 GDP
## [1,] 9.378173e-01 1.930019e-01 -2145.1937095
## [2,] 1.349227e+00 -3.012665e-02 -437.2668661
## [3,] -1.527117e-01 4.467480e-01 738.9138900
## [4,] 4.599431e-06 2.622808e-05 0.6212176
## [5,] -4.605570e-01 2.031003e-02 694.1695746
## [6,] 3.501744e-01 -5.858274e-02 -446.3634874
## [7,] -7.664498e-06 -2.396556e-05 0.4071503print(Sigma_hat)## UEq CPIq GDP
## UEq 2.528039e-01 -0.006658128 -161.21120
## CPIq -6.658128e-03 0.068874255 -71.92108
## GDP -1.612112e+02 -71.921082111 3579281.40727Cholesky Factorization in R
Sigma_hat_tr <- t(chol(Sigma_hat))
print(Sigma_hat_tr)## UEq CPIq GDP
## UEq 0.5027961 0.0000000 0.000
## CPIq -0.0132422 0.2621047 0.000
## GDP -320.6294114 -290.5972613 1841.747# Verify that Sigma_hat_tr*t(Sigma_hat_tr) = Sigma_hat
print(Sigma_hat)## UEq CPIq GDP
## UEq 2.528039e-01 -0.006658128 -161.21120
## CPIq -6.658128e-03 0.068874255 -71.92108
## GDP -1.612112e+02 -71.921082111 3579281.40727print(Sigma_hat_tr %*% t(Sigma_hat_tr))## UEq CPIq GDP
## UEq 2.528039e-01 -0.006658128 -161.21120
## CPIq -6.658128e-03 0.068874255 -71.92108
## GDP -1.612112e+02 -71.921082111 3579281.40727Procedure 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="GDP Response", col="blue", lwd=2)
abline(h=0,col="black",lwd=1)