In [1]:
%pylab inline
This worksheet demonstrates how to solve the RBC model with the dolo library and how to exploit the results.
examples\notebooks\. The notebook was opened and run from that directory.examples\global_models\
In [2]:
# First we import the dolo library
from dolo import *
In [3]:
# The rbc model is defined in a yaml file (see doc for syntax)
!cat '../global_models/rbc.yaml'
In [4]:
model = yaml_import('../global_models/rbc.yaml')
In [5]:
# let have a look at the model equations and check that
# all residual errors are 0 at the steady-state
print( model ) # sympified output is currently broken
In [6]:
dr_pert = approximate_controls(model, order=1)
#dr = approximate_controls(model, order=2)
In [7]:
dr_global = time_iteration(model, pert_order=1, )
In [8]:
bounds = [dr_global.smin[1], dr_global.smax[1]]
plot_decision_rule(model, dr_global, 'k', 'i', label='global', bounds=bounds)
plot_decision_rule(model, dr_pert, 'k', 'i', label='perturbation', bounds=bounds)
ylabel('i')
title('Investment')
legend()
show()
display("Maybe the approximation method doesn't make a big difference in this model...")
Repeat the same exercise for various parameter values
In [9]:
model.symbols['states']
Out[9]:
In [10]:
model.calibration['parameters'] # let check current parameter values
Out[10]:
In [11]:
drs = []
values = linspace(0.01, 0.04,5)
for val in values:
model.set_calibration(delta=val)
drs.append(approximate_controls(model))
In [12]:
for i,dr in enumerate(drs):
plot_decision_rule(model, dr, 'k', 'i', label='$\delta={}$'.format(values[i]), bounds=bounds)
ylabel('i')
title('Investment')
legend()
show()
In [13]:
# we will use the deterministic steady-state as a starting point
s0 = model.calibration['states']
# we also get the covariance matrix just in case
sigma = model.get_distribution().sigma
In [14]:
s1 = s0.copy()
s1[0] += 0.1
In [15]:
irf = simulate(model, dr_global, s1, sigma, n_exp=0, horizon=40 )
# choosing n_exp>1, will result in many "stochastic" simulations
# with n_exp = 0, we get one single simulation without shock
# the output is a panda table
print(irf.__class__)
print(irf.shape)
# the result is of size n_v x H where n_v is the number of variables
# in the model and H the number of dates
In [16]:
# let plot the response of consumption and investment
figsize(10,10)
subplot(221)
plot(irf['z'], label='productivity')
subplot(222)
plot(irf['i'], label='investment')
subplot(223)
plot(irf['n'], label='labour')
subplot(224)
plot(irf['c'], label='consumption')
Out[16]:
Note that the plotting is made using the wonderful matplotlib library. See the doc online to customize the plots to your needs.
In [17]:
# export to matlab
# we need to convert the dataframe to an array:
irf_array = array( irf )
import scipy.io
scipy.io.savemat("export.mat", {'table': irf_array} )
In [18]:
sim = simulate(model, dr_global, s0, sigma, n_exp=1000, horizon=40 )
# now we run 1000 random simulations
# the result is an array of size H x n_v x n_exp where
# - H the number of dates
# - n_v is the number of variables
# - n_exp the number of simulations
print(sim.shape)
In [19]:
sim.shape
Out[19]:
In [20]:
# let plot the response of consumption and investment
i_z = model.variables.index('z')
i_i = model.variables.index('i')
i_n = model.variables.index('n')
i_c = model.variables.index('c')
figsize(10,10)
for i in range(1000):
subplot(221)
plot(sim[:, i, i_z], color='red', alpha=0.1)
subplot(222)
plot(sim[:, i, i_i], color='red', alpha=0.1)
subplot(223)
plot(sim[:, i, i_n], color='red', alpha=0.1)
subplot(224)
plot(sim[:, i, i_c], color='red', alpha=0.1)
subplot(221)
title('productivity')
subplot(222)
title('investment')
subplot(223)
title('labour')
subplot(224)
title('consumption')
Out[20]:
In [21]:
# We can also make some statistics
dsim = log(sim[1:,:,:]/sim[:-1,:,:,]) # compute growth rates
volat = dsim.std(axis=0) # compute volatility of growth rates for each simulation
table = column_stack([
volat.mean(axis=0),
volat.std(axis=0)
])
display(table)
print('First column contains the standard deviations of growth rates.')
print('Second column contains the confidence intervals.')
We can use the pandas library to present the results
In [22]:
model.variables
Out[22]:
In [24]:
import pandas
df = pandas.DataFrame(table, index=model.variables[:-1], columns=['std. deviations','confidence intervals'])
df
Out[24]:
In [25]:
# export to matlab
# - stochastic simulations are already stored as an array
# - the dataframe needs to be converted
import scipy.io
scipy.io.savemat("export2.mat", {'stoch_sim': sim, 'stats': array(df)} )
In [26]:
from dolo.numeric.error_measures import omega
In [27]:
bounds = row_stack( [dr_global.smin, dr_global.smax] )
display(bounds)
In [29]:
omega(dr, model, bounds, [50,50], time_weight=[100, 0.96, dr_pert.S_bar])
Out[29]:
In [ ]: