Loading packages etc.

We first clear the workspace using rm(list = ls())and then include all packages we need. If a package is missing in your R distribution (which is quite likely initially), just use install.packages("package_name") with the respective package name to install it on your system. If you execute the code in the file install_packages.R, then all necessary packages will be installed into your R distribution. If the variable export_graphs is set to TRUE, then the graphs will be exported as pdf-files.

rm(list = ls())
library(reshape2)
library(base)
library(ggplot2)
library(grid)
library(stringr)
library(tidyverse)

# should graphs be exported to pdf
export_pdf = FALSE

Plot a time path of consumption

Before plotting a time path of consumption, we have to parameterize the problem, i.e. we have to choose \(\gamma\), \(\beta\) and \(a_0\). We set some choices here that can easily be changed.

# curvature of utility function
gamma = 0.5

# time discount factor
beta = 0.95

# initial quantity of the resource
a0 = 100


Now we can calculate a time path using the explicit solution we obtained. We therefore choose to simulate 200 periods. In the data frame, we store the different period dates as well as the corresponding consumption level \(\beta^{t\gamma}\cdot(1-\beta^\gamma)\cdot a_0\). We then plot the consumption path using ggplot.

# simulate path for 200 periods
t = c(0:200)

# generate data frame for plotting
dat <- data.frame(
  t = t,
  c = beta^(t*gamma)*(1-beta^gamma)*a0
)

# now generate time path plot
myplot <- ggplot(data=dat) + 
  geom_line(aes(x=t, y=c), color="darkblue", linewidth=1) +
  scale_x_continuous(breaks=seq(0, 200, 20)) +
  labs(x = expression(paste("Time ", t)),
       y = expression(paste("Consumption path ", c[t]))) +
  theme_bw()

# print plot
print(myplot)

Value and policy function

Next, we would like to draw the policy and the value function of the cake eating problem. We therefore define a discrete set of points \(a \in [0, a_0]\) for which we calculate the policy and value function values. We can do this using the explicit solutions provided in class. Finally, we can plot the policy function by means of ggplot.

# generate discrete set for potential resources
a = seq(0, a0, 0.01)

# generate data frame for plotting
dat <- data.frame(
  a = a,
  c = a*(1-beta^gamma),
  V = (1-beta^gamma)^(-1/gamma)*a^(1-1/gamma)/(1-1/gamma)
)

# plot for the policy function
myplot <- ggplot(data=dat) + 
  geom_line(aes(x=a, y=c), color="darkblue", linewidth=1) +
  scale_x_continuous(breaks=seq(0, a0, 20)) +
  labs(x = expression(paste("Resource ", a[t])),
       y = expression(paste("Policy function ", c(a[t])))) +
  theme_bw()

# print plot
print(myplot)


And we can do the same for the value function.

# plot for the value function
myplot <- ggplot(data=dat) + 
  geom_line(aes(x=a, y=V), color="darkred", linewidth=1) +
  coord_cartesian(ylim=c(-1600, 0)) + 
  scale_x_continuous(breaks=seq(0, a0, 20)) +
  labs(x = expression(paste("Resource ", a[t])),
       y = expression(paste("Value function ", V(a[t])))) +
  theme_bw()

# print plot
print(myplot)

Simulating the system forward

Now we want to use the policy function \(c(a_t)\) as well as the low of motion for the resource \(a_{t+1} = a_t - c_t\) to simulate the solution of the dynamic program forward. To do so, we first define a function policy that returns the value of the policy function for any given state a. Next to the current level of the resource, this function receives the parameters of the problem gamma and beta.

policy <- function(a, gamma, beta) {
  
  # calculate the optimal policy
  c = (1-beta^gamma)*a
  
  # return the value
  return(c)
}


Now we can sequentially simulate the system forward. We start by defining a final period T until which we want to simulate the system. We set the inital state of the resource to a_t = a0 and calculate the optimal policy at this point. We then successively use the budget constraint and the optimal policy function to calculate next period’s resource level and next period’s optimal consumption until we reach the final simulation period T. Finally, we can plot the optimal consumption path as well as the optimal path as well as the optimal path for the evolution of the resource over time.

# number of periods to simulate
T = 200

# start with the initial value of the resource
a_t = a0

# calculate the policy at this point
c_t = policy(a_t[1], gamma, beta)

# now simulate forward until date T
for(t in 2:(T+1)) {
  a_t[t] = a_t[t-1] - c_t[t-1]
  c_t[t] = policy(a_t[t], gamma, beta)
}

# the time periods
t = c(0:200)

# generate data frame for plotting
dat <- data.frame(
  t = t,
  c = c_t,
  a = a_t
)

# now generate time path plot of consumption
myplot <- ggplot(data=dat) + 
  geom_line(aes(x=t, y=c), color="darkblue", linewidth=1) +
  scale_x_continuous(breaks=seq(0, 200, 20)) +
  labs(x = expression(paste("Time ", t)),
       y = expression(paste("Consumption path ", c[t]))) +
  theme_bw()

print(myplot)

# now generate time path plot of the resource
myplot <- ggplot(data=dat) + 
  geom_line(aes(x=t, y=a), color="darkred", linewidth=1) +
  scale_x_continuous(breaks=seq(0, 200, 20)) +
  labs(x = expression(paste("Time ", t)),
       y = expression(paste("Resource level ", a[t]))) +
  theme_bw()

print(myplot)