Computing Economic Models in R: How to apply shocks to parameter values in the euler equation?

Question

Hi everyone im using R to try and simulate some economic models. We do this primarily through the use of the euler equation. I've figured out that applying shocks to values which are defined within the function (in this case it is k is pretty simple as seen in the code below, however I'm interested in applying a shock to parameters like delta, theta and rho.

For what its worth I'm using the R package deSolve. Any help is appreciated.

library('deSolve')

##############################################
#Computing the neoclassical growth model in R#
##############################################

#parameters and state space
A<-1
theta<- 0.1
alpha<-0.5
delta<-0.3
rho<-0.9
kinital <- c(k = 1)
times <- seq(from = 0, to = 100, by = 0.2)

#define euler equation
euler <- function(t, k, parms)
  list((1/theta)*alpha*A*k^(alpha-1)-delta-rho)

#Compute 
out <- ode(y = kinital, times = times, func = euler,
           parms = NULL)

plot(out, main = "Euler equation", lwd = 2)


#########################
#Temporary Capital Shock#
########################
eventdat <- data.frame(var = c("k"),
                       time = c(30) ,
                       value = c(10),
                       method = c("add"))
eventdat1 <- data.frame(var = c("k"),
                       time = c(30) ,
                       value = c(-5),
                       method = c("add"))
out3<-ode(y=kinital,times=times,func=euler,events=list(data=eventdat))
out4<-ode(y=kinital,times=times,func=euler,events=list(data=eventdat1))

plot(out,out3,out4,main="Temporary Shock",lwd=3)

Answer 1

A colleague off of stackexchange suggested a cleaner bit of code which is a bit cleaner. This is seen below:

A<-1
theta<- 0.1
alpha <- 0.5
rho<-0.9
init <- c(k = 17, delta = 0.3)
times <- seq(from = 0, to = 400, by = 0.2)

euler.function<-function(t,y, prams){
  k <- y[1]
  delta <- y[2]
  dk <- (1/theta)*alpha*A*k^(alpha-1)-delta-rho
  list(c(dk, 0))}

deventdat<- data.frame(var = c("delta", "delta"),
                       time = c(30, 51) ,
                       value = c(0.1, -0.1),
                       method = c("add"))

res<-ode(y=init,times=times, func=euler.function, parms=NULL, events=list(data=deventdat))

plot(res,lwd=2)

Answer 2

Not a great fix but the way to deal with this type of problem is by conditioning your values to take place over some interval. I do this for depreciation as follows:

##############################
#Temporary Depreciation Shock#
##############################
#New Vars
A<-1
theta<- 0.1
alpha<-0.5
delta<-0.3
rho<-0.9
kinital <- c(k = 17)
times <- seq(from = 0, to = 400, by = 0.2)
#Redefine Euler
euler2<-function(t,k,prams){
  list((1/theta)*alpha*A*k^(alpha-1)-delta-rho)}

euler3<-function(t,k,prams){
  list((1/theta)*alpha*A*k^(alpha-1)-(delta+0.05*(t>=30&t<=40))-rho)}


#Output
doutbase<-ode(y=kinital,times=times, func=euler2, parms=NULL)
doutchange<-ode(y=kinital,times=times, func=euler3, parms=NULL)

#plots
plot(doutbase,doutchange,main="Change in depreciation at t=30 until t=40",lwd=2)