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)

You can use a forcing function instead of an event, see [example here](https://tpetzoldt.github.io/deSolve-forcing/deSolve-forcing.html). More examples can be found on StackOverflow, just search for [`[r]` desolve forcing](https://stackoverflow.com/search?q=%5Br%5D+desolve+forcing). Another method is to use a "dummy" state variable with derivative zero and then apply events. — tpetzoldt, May 05 '22 at 19:39

EconJohn · Answer 1 · 2022-05-06T02:45:23.720

1

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)

edited May 06 '22 at 02:45

answered May 06 '22 at 01:04

EconJohn

165
11

Such "if" or "<=" approaches are simple and work in practice, but are not ideal for a smooth differential equation system. – tpetzoldt May 06 '22 at 15:55
@tpetzoldt can you elaborate or direct me to a source for why not? – EconJohn May 06 '22 at 15:57
We had this discussion already on SO (and a past useR! conference), but it can be explained with simple words: An ODE is **continuous** and differentiable by definition. A discontinuity "irritates" the solvers, so that they may iterate a lot in front of it, like a car hitting a kerbstone. The good news is, that deSolve's solvers are quite robust in practice. In effect, discontinuities cost CPU time. Events stop the solver at the discontinuity, change the states and then restart. This is cleaner, but costs also CPU time and is programmatically more complex. The final answer is: "it depends". – tpetzoldt May 06 '22 at 16:05
@tpetzoldt thank you so much for the detailed answer! I havent run into problems on my machine so far, but as I expand these models Im sure CPU time will get expensive. Thank you! – EconJohn May 06 '22 at 18:17

EconJohn · Accepted Answer · 2022-05-06T04:20:54.077

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)

edited May 06 '22 at 04:20

answered May 06 '22 at 03:45

EconJohn

165
11

1

Yes, this is essentially the approach with an additional "dummy" state variable that is controlled by an event. More complex than the "<=" approach and IMHO "philosophically" better. – tpetzoldt May 06 '22 at 15:58

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

2 Answers2