Final Paper: Testing Code on Simulated Data

Set work directory

setwd("C:/Users/victo/Downloads/ECON 4310/Exercise")

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.3

US Data From Class

data    <- read.csv("R_UR_PCE.csv", header=TRUE)
FFR     <- ts(data$FEDFUNDS, frequency=12, start=c(1959,1))
UR      <- ts(data$UNRATE, frequency=12, start=c(1959,1))
PCE     <- ts(data$PCEPILFE, frequency=12, start=c(1959,1))
INFL    <- 1200*diff(log(PCE))
FFR     <- window(FFR,c(1960,1),c(2020,1))
UR      <- window(UR,c(1960,1),c(2020,1))
INFL    <- window(INFL,c(1960,1),c(2020,1))

Estimate VAR

# Prepare data
FFR  <- window(FFR,c(1984,1),c(2006,4))
UR   <- window(UR,c(1984,1),c(2006,4))
INFL <- window(INFL,c(1984,1),c(2006,4))
Tall <- length(INFL) # total number of observations, starting 1984:M1
T1   <- (1+12) # estimation sample will start in T1 = 1985:Q1
p    <- 2 # number of lags
n    <- 3
yyall <- cbind(FFR,UR,INFL)

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),]
  
}

Compute 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)

##              FFR           UR        INFL
## [1,]  0.15078450  0.108953815 -1.70207812
## [2,]  1.34436825 -0.201572713  0.55127471
## [3,] -0.32194100  0.791665697  0.52755007
## [4,]  0.02134402  0.004869783  0.01666898
## [5,] -0.36252381  0.199668669 -0.26254317
## [6,]  0.29417067  0.181697060 -0.03099540
## [7,]  0.01442170  0.009764637 -0.04721704

print(Sigma_hat)

##               FFR           UR         INFL
## FFR   0.034323148 -0.006268862  0.037896164
## UR   -0.006268862  0.016746806 -0.001034543
## INFL  0.037896164 -0.001034543  1.824616397

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)){
    epst <- rnorm(n, mean=1, sd=1)
    ut <- Sigma_tr*epst
    yy[t,] <- xt%*%Phi + ut
    # 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)
  }

Simulations for Forecasting

ForecastVARp <- 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)){
    epst <- rnorm(n, mean=1, sd=1)
    ut <- Sigma_tr*epst
    yy[t,] <- xt%*%Phi + ut
    # 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)
}

# recursive estimation
# forecast generation