# Example data generation: Y is an AR(1) process that may receive a shock at
# t=50, S is the shock (0/1), a combination of 3 AR(1) processes (X1-X3)
# X4 is another AR(1) process, uncorrelated with S, X4sq is just X4^2
# All AR(1) processes have the same phi=0.98 coefficient, and are Monte
# Carlo sampled 500 times
set.seed( 1 ) # make results reproducible
# LSD-like arrays to store simulated time series (t x var x MC)
dataNoShock <- dataShock <-array ( 0, dim = c( 60, 7, 500 ) )
colnames( dataNoShock ) <- colnames( dataShock ) <-
c( "Y", "S", "X1", "X2", "X3", "X4", "X4sq" )
# Monte Carlo sampling
for( n in 1 : 500 ) {
# simulation time
for( t in 2 : 60 ) {
# AR process on X vars
for( v in c( "X1", "X2", "X3", "X4" ) ) {
dataNoShock[ t, v, n ] = dataShock[ t, v, n ] =
0.98 * dataShock[ t - 1, v, n ] + rnorm( 1, 0, 0.1 )
}
# apply shock once
if( t == 50 ) {
dataShock[ t, "S", n ] <- 1
shockEff <- 0.4 + 0.7 * isTRUE( dataShock[ t, "X1", n ] > 0.1 ) -
0.4 * isTRUE( dataShock[ t, "X2", n ] > 0.1 ) +
0.2 * isTRUE( dataShock[ t, "X3", n ] > 0.05 ) + rnorm( 1, 0, 0.2 )
} else
shockEff <- 0
# AR process on Y var
rs <- rnorm( 1, 0, 0.1 )
dataNoShock[ t, "Y", n ] = 0.98 * dataNoShock[ t - 1, "Y", n ] + rs
dataShock[ t, "Y", n ] = 0.98 * dataShock[ t - 1, "Y", n ] + shockEff + rs
}
}
# another uncorrelated var
dataNoShock[ , "X4sq", ] <- dataShock[ , "X4sq", ] <- dataShock[ , "X4", ] ^ 2
# \donttest{
# linear IRF analysis
linearIRF <- irf.lsd( data = dataNoShock, # non-shocked MC data
data.shock = dataShock, # shocked data
t.horiz = 10, # post-shock analysis time horizon
var.irf = "Y", # variable to compute IRF
var.shock = "S" ) # shock variable (impulse)
plot( linearIRF, irf.type = "cum.irf" ) # cumulative IRF plot
print( linearIRF ) # show IRF data
# }
Run the code above in your browser using DataLab