Here is my code:
clear
clc
close all
toolkit_folder = 'C:\Users\aledi\Documents\GitHub\VFIToolkit-matlab';
addpath(genpath(toolkit_folder))
%% Begin setting up to use VFI Toolkit to solve
% Lets model agents from age 20 to age 100, so 81 periods
Params.agejshifter=19; % Age 20 minus one. Makes keeping track of actual age easy in terms of model age
Params.J=100-Params.agejshifter; % =81, Number of period in life-cycle
% Grid sizes to use
n_d = 0; % Endogenous labour choice (fraction of time worked)
d_grid = 0;
n_a = 301; % Endogenous asset holdings
n_z = [7,5]; % Exogenous labor productivity units shocks, persistent and transitiory
N_j = Params.J; % Number of periods in finite horizon
%% Parameters
% Discount rate
Params.beta = 0.96;
% Preferences
Params.crra = 2; % Coeff of relative risk aversion (curvature of consumption)
% Prices
Params.w=1; % Wage
Params.r=0.05; % Interest rate (0.05 is 5%)
% Demographics
Params.agej=1:1:Params.J; % Is a vector of all the agej: 1,2,3,...,J
Params.Jr=46;
% Pensions
Params.pension=0.3;
% Age-dependent labor productivity units
Params.kappa_j=[linspace(0.5,2,Params.Jr-15),linspace(2,1,14),zeros(1,Params.J-Params.Jr+1)];
% persistent AR(1) process on idiosyncratic labor productivity units
rho_z1=0.9;
sigma_epsilon_z1=0.02;
% transitory iid normal process on idiosyncratic labor productivity units
sigma_epsilon_z2=0.2; % Implictly, rho_z2=0
% Conditional survival probabilities: sj is the probability of surviving to be age j+1, given alive at age j
% Most countries have calculations of these (as they are used by the government departments that oversee pensions)
% In fact I will here get data on the conditional death probabilities, and then survival is just 1-death.
% Here I just use them for the US, taken from "National Vital Statistics Report, volume 58, number 10, March 2010."
% I took them from first column (qx) of Table 1 (Total Population)
% Conditional death probabilities
Params.dj=[0.006879, 0.000463, 0.000307, 0.000220, 0.000184, 0.000172, 0.000160, 0.000149, 0.000133, 0.000114, 0.000100, 0.000105, 0.000143, 0.000221, 0.000329, 0.000449, 0.000563, 0.000667, 0.000753, 0.000823,...
0.000894, 0.000962, 0.001005, 0.001016, 0.001003, 0.000983, 0.000967, 0.000960, 0.000970, 0.000994, 0.001027, 0.001065, 0.001115, 0.001154, 0.001209, 0.001271, 0.001351, 0.001460, 0.001603, 0.001769, 0.001943, 0.002120, 0.002311, 0.002520, 0.002747, 0.002989, 0.003242, 0.003512, 0.003803, 0.004118, 0.004464, 0.004837, 0.005217, 0.005591, 0.005963, 0.006346, 0.006768, 0.007261, 0.007866, 0.008596, 0.009473, 0.010450, 0.011456, 0.012407, 0.013320, 0.014299, 0.015323,...
0.016558, 0.018029, 0.019723, 0.021607, 0.023723, 0.026143, 0.028892, 0.031988, 0.035476, 0.039238, 0.043382, 0.047941, 0.052953, 0.058457, 0.064494,...
0.071107, 0.078342, 0.086244, 0.094861, 0.104242, 0.114432, 0.125479, 0.137427, 0.150317, 0.164187, 0.179066, 0.194979, 0.211941, 0.229957, 0.249020, 0.269112, 0.290198, 0.312231, 1.000000];
% dj covers Ages 0 to 100
Params.sj=1-Params.dj(21:101); % Conditional survival probabilities
Params.sj(end)=0; % In the present model the last period (j=J) value of sj is actually irrelevant
figure
plot(Params.agej,Params.kappa_j)
hold on
plot(Params.agej,Params.pension*ones(N_j,1))
%% Grids
% The ^3 means that there are more points near 0 and near 10. We know from
% theory that the value function will be more 'curved' near zero assets,
% and putting more points near curvature (where the derivative changes the most) increases accuracy of results.
a_min = 0;
a_max = 10;
a_grid=a_min+(a_max-a_min)*(linspace(0,1,n_a).^3)'; % The ^3 means most points are near zero, which is where the derivative of the value fn changes most.
% First, the AR(1) process z1
Tauchen_q = 3.0;
[z1_grid,pi_z1]=discretizeAR1_Tauchen(0.0,rho_z1,sigma_epsilon_z1,n_z(1),Tauchen_q);
z1_grid = exp(z1_grid);
% Second, the transitory process. Since z2 is purely iid, all rows of the
% transition matrix pi_z2 are the same, i.e. Prob(z2,z2') does not depend
% on z2.
[z2_grid,pi_z2]=discretizeAR1_Tauchen(0.0,0.0,sigma_epsilon_z2,n_z(2),Tauchen_q);
z2_grid = exp(z2_grid);
pi_e = pi_z2(1,:)'; % Distribution of z2=e, transitory shock
mean_z2 = dot(pi_e,z2_grid);
z2_grid = z2_grid/mean_z2;
% Combine transition matrix for two shocks. z1 is ordered first, Kronecker
% in reverse order.
z_grid = [z1_grid;z2_grid];
pi_z = kron(pi_z2,pi_z1);
%% Set up return function
DiscountFactorParamNames={'beta','sj'};
% input z1 and z2
ReturnFn=@(aprime,a,z1,z2,agej,Jr,w,kappa_j,r,pension,crra) ...
Model_ReturnFn(aprime,a,z1,z2,agej,Jr,w,kappa_j,r,pension,crra);
%% Now solve the value function iteration problem, just to check that things are working before we go to General Equilbrium
disp('Test ValueFnIter')
vfoptions=struct(); % Just using the defaults.
vfoptions.verbose = 1;
vfoptions.divideandconquer = 0;
tic;
[V, Policy]=ValueFnIter_Case1_FHorz(n_d,n_a,n_z,N_j,d_grid,a_grid,z_grid,pi_z,ReturnFn,Params,DiscountFactorParamNames,[],vfoptions);
toc
% V is now (a,z1,z2,j). One dimension for each state variable.
% Compare
size(V)
% with
[n_a,n_z(1),n_z(2),N_j]
% there are the same.
% Policy is
size(Policy)
% which is the same as
[length(n_a),n_a,n_z(1),n_z(2),N_j]
% The n_a,n_z(1),n_z(2),N_j represent the state on which the decisions/policys
% depend, and there is one decision for each decision variable 'd' and each
% endogenous state variable 'a', and one for each exogenous state variable 'z'
If I set vfoptions.divideandconquer = 0, or if I leave the default options, all good. But if I set vfoptions.divideandconquer = 1, then I get the error:
Error using CreateReturnFnMatrix_Case1_Disc_DC1_nod_Par2
Too many input arguments.
Error in ValueFnIter_Case1_FHorz_DC1_nod_raw (line 24)
ReturnMatrix_ii=CreateReturnFnMatrix_Case1_Disc_DC1_nod_Par2(ReturnFn, n_z, a_grid, a_grid(level1ii), z_gridvals_J(:,:,N_j), ReturnFnParamsVec,1);
