Solving dynamic programming models with the NEGM

This notebooks shows how the model from the introduction of "A Fast Nested Endogenous Grid Method for Solving General Consumption-Saving Models" (Druedahl, 2017) is solved in just a few seconds. Basic knowledge og MATLAB and C++ is presumed.
(Code for two-asset model with bells and whistles http://www.econ.ku.dk/druedahl/code/NEGM/code_twoasset.zip):
This code package requires:
  1. Windows 64-bit
  2. MATLAB 2017b (not tested with earlier versions).
  3. GCC 5.3.0 (installation covered below)
  4. NLopt 2.4.2 C++ optimizer (https://nlopt.readthedocs.io/, normally no installation required)
All comments and questions should be adressed to jeppe.druedahl@econ.ku.dk.
% first let us clear all memory and close all figures
clc;
clear;
close all;

0. Model

To recap the model is:
model.png
The model is solved in ratio form. Lowercase variables denote normalization with permanent income. Households live for T periods. In the first period they work and . In the remaining periods they are retired and and . Details are in the paper.

1. Installation of GCC (and NLopt)

To compile C++ files it is necessary to install a compiler. The GCC compiler is free and can be installed by:
  1. In the Apps ribbon choose "Get More Apps"
  2. Locate "MATLAB Support of MinGW-w64 C/C++ Compiler"
  3. Install following the instructions
Ensure that MinGW64 Compiler (C++) is used:
% run and choose MinGW64 if not already chosen
mex -setup C++
MEX configured to use 'MinGW64 Compiler (C++)' for C++ language compilation. Warning: The MATLAB C and Fortran API has changed to support MATLAB variables with more than 2^32-1 elements. You will be required to update your code to utilize the new API. You can find more information about this at: http://www.mathworks.com/help/matlab/matlab_external/upgrading-mex-files-to-use-64-bit-api.html. To choose a different C++ compiler, select one from the following: MinGW64 Compiler (C++) mex -setup:C:\Users\okoJD\AppData\Roaming\MathWorks\MATLAB\R2017b\mex_C++_win64.xml C++ Microsoft Visual C++ 2015 mex -setup:'C:\Program Files\MATLAB\R2017b\bin\win64\mexopts\msvcpp2015.xml' C++
SKIP AT FIRST: If you experience errors related to NLopt it probably needs to be recompilled for your system as follows:
  1. Install TDM-GCC from http://tdm-gcc.tdragon.net/
  2. Open "MinGW Command Prompt" and relocate to "cfuncs/nlopt-2.4.2-dll64/"
  3. Run "dlltool --input-def libnlopt-0.def --dllname libnlopt-0.dll --output-lib libnlopt-0.lib"
  4. Copy "libnlopt-0.dll" to the current folder

2. Overview

Let us get a short overview of the code base.
Current folder:
cfuncs/: all the C++ functions
figs/: stores figure output

2. Compile mex

To compile the C++ files so that they can be called from MATLAB the following lines are necessary:
threads = 8; % max number of threads used in parallization
model.mex_solve(threads)
Building with 'MinGW64 Compiler (C++)'. MEX completed successfully.
model.mex_simulate(threads)
Building with 'MinGW64 Compiler (C++)'. MEX completed successfully.

3. Setup, Solve and Simulate the model

% a. setup baseline parameters
par = model.setup();
% b. solve the model
t0 = tic;
[par,sol] = model.solve(par);
fprintf('time to solve: %5.2f secs',toc(t0));
time to solve: 1.83 secs
% c. simulate the model
t0 = tic;
sim = model.simulate(par,sol);
fprintf('time to simulate: %5.2f secs',toc(t0));
time to simulate: 0.65 secs

4. Figures

funs.layout();
t = 1;
Nb_vec = floor([par.Ndb/5 par.Ndb/4 par.Ndb/3 par.Ndb/2]);
% a. consumption functions
fig = figure('name','consumption_function');
hold on;
for i_db = Nb_vec
plot(par.grid_m+par.grid_db(i_db),sol.c{t,1}(:,i_db),'-',...
'DisplayName',sprintf('keep, $\\bar{d}_t = %5.2f$',par.grid_db(i_db)));
end
plot(par.grid_x,sol.c{t,2},'-','color','black','DisplayName','adj.');
grid on;
ylabel('$c_t$')
xlabel('$x_t = m_t + \bar{d}_t$')
xlim([0 5])
legend('show','Location','best')
funs.printfig(fig);
% b. durable and saving adjuster function
fig = figure('name','d_a_adj_function');
hold on;
plot(par.grid_x,sol.d{t,2},'-','DisplayName','$d_t$');
plot(par.grid_x,par.grid_x-sol.d{t,2}-sol.c{t,2},'-','DisplayName','$a_t$');
grid on;
ylabel('$d_t$')
xlabel('$x_t = m_t + \bar{d}_t$')
xlim([0 5])
legend('show','Location','best')
funs.printfig(fig);
% c. value functions
fig = figure('name','value_function');
hold on;
for i_db = Nb_vec
plot(par.grid_m+par.grid_db(i_db),sol.v{t,1}(:,i_db),'-',...
'DisplayName',sprintf('keep, $\\bar{d}_t = %5.2f$',par.grid_db(i_db)));
end
plot(par.grid_x,sol.v{t,2},'color','black','DisplayName','adj.');
grid on;
ylabel('$v_t$')
xlabel('$x_t = m_t + \bar{d}_t$')
xlim([0 5])
legend('show','Location','best')
funs.printfig(fig)
% d. life-cycle profiles
fig = figure('name','life_cycle_profiles');
hold on;
plot(par.age(10:par.TR),mean(sim.c(:,10:par.TR)),'-','DisplayName','$c_t$');
plot(par.age(10:par.TR),mean(sim.d(:,10:par.TR)),'-','DisplayName','$d_t$');
plot(par.age(10:par.TR),mean(sim.a(:,10:par.TR)),'-','DisplayName','$a_t$');
grid on;
xlabel('age')
legend('show','Location','best')
funs.printfig(fig);
% e. share of adjusters
fig = figure('name','share_of_adjusters');
hold on;
plot(par.age(10:par.TR),mean(sim.z(:,10:par.TR)),'-');
grid on;
xlabel('age')
funs.printfig(fig);
% f. histograms
fig = figure('name','histogram_m_age_35_65');
histogram(funs.vec(sim.m(:,10:par.TR)),100,'normalization','probability');
xlabel('$m_t$')
grid on;
funs.printfig(fig);
fig = figure('name','histogram_db_age_35_65');
histogram(funs.vec(sim.db(:,10:par.TR)),100,'normalization','probability');
xlabel('$\bar{d}_t$')
grid on;
funs.printfig(fig);