[This post is a continuation of Part I]
In the previous post I outline a solution that requires to change the code in CreateReturnFnMatrix_Case1_Disc and change Aiyagari1994_ReturnFn as well. This requires to keep two versions of Aiyagari1994_ReturnFn, one for the GPU (not vectorized) and one for the CPU (vectorized). To avoid this, I tried to use the option vfoptions.returnmatrix==1 suggested by the Toolkit.
As far as I understand, I have to create myself the 3-dimensional return matrix over (a’,a,z) and pass it to VFI as the input ReturnFn. In the standard case, ReturnFn is a function handle, but now it is a numeric array (right?). So here’s my code:
% Pre-compute here the return matrix
tic
if vfoptions.parallel==0 && vfoptions.returnmatrix==1
z_grid = gather(z_grid);
pi_z = gather(pi_z);
ReturnFn = zeros(n_a,n_a,n_z);
for i3=1:n_z
for i2=1:n_a
for i1=1:n_a
ReturnFn(i1,i2,i3)=Aiyagari1994_ReturnFn(a_grid(i1),a_grid(i2),z_grid(i3),Params.alpha,Params.delta,Params.mu,Params.r);
end
end
end
end
toc
This is quite fast even if it is not vectorized. It takes 0.6 seconds with n_a=512, n_z=11.
Unfortunately, I got the following error message
Incorrect number or types of inputs or outputs for function func2str.
Error in getAnonymousFnInputNames (line 3)
temp=func2str(FnToEvaluate); % Note that this strips out any spaces automatically
Error in ValueFnIter_Case1 (line 191)
temp=getAnonymousFnInputNames(ReturnFn);
Error in Aiyagari1994 (line 130)
[V,Policy]=ValueFnIter_Case1(n_d,n_a,n_z,d_grid,a_grid,z_grid,pi_z,ReturnFn,Params,DiscountFactorParamNames,[],vfoptions);
For completeness, I report my entire code below:
% Example based on Aiyagari (1994).
%
% These codes set up and solve the Aiyagari (1994) model for a given
% parametrization. After solving the model they then show how some of the
% vfitoolkit commands to easily calculate things like the Gini coefficient
% for income, and how to plot the distribution of asset holdings.
%
% VFI Toolkit automatically detects hardware (GPU? Number of CPUs?) and
% sets defaults accordingly. It will run without a GPU, but slowly. It is
% indended for use with GPU.
clear,clc,close all
addpath(genpath('C:\Users\aledi\Dropbox\Family Room\6_Computation\VFI_toolkit\VFIToolkit-matlab-master'))
%% Set some basic variables
% VFI Toolkit thinks of there as being:
% k: an endogenous state variable (assets)
% z: an exogenous state variable (exogenous labor supply)
% Size of the grids
n_k = 512;
n_z = 11; %21;
% Parameters
Params.beta = 0.96; %Model period is one-sixth of a year
Params.alpha = 0.36;
Params.delta = 0.08;
Params.mu = 3;
Params.sigma = 0.2;
Params.rho = 0.6;
%% Set up the exogenous shock process
% Create markov process for the exogenous labour productivity, l.
Tauchen_q=3; % Footnote 33 of Aiyagari(1993WP, pg 25) implicitly says that he uses q=3
[z_grid,pi_z]=discretizeAR1_Tauchen(0,Params.rho,sqrt((1-Params.rho^2)*Params.sigma^2),n_z,Tauchen_q);
% Note: sigma is standard deviations of s, input needs to be standard deviation of the innovations
% Because s is AR(1), the variance of the innovations is (1-rho^2)*sigma^2
z_grid = exp(z_grid);
aux = pi_z^1000;
z_prob = aux(1,:)';
z_mean = dot(z_grid,z_prob);
z_grid = z_grid/z_mean;
%[z_mean,z_variance,z_corr,~]=MarkovChainMoments(z_grid,pi_z);
% Get some info on the markov process
%[Expectation_l,~,~,~]=MarkovChainMoments(z_grid,pi_z); %Since l is exogenous, this will be it's eqm value
% Note: Aiyagari (1994) actually then normalizes l by dividing it by Expectation_l (so that the resulting process has expectation equal to 1)
%z_grid=z_grid./Expectation_l;
%[Expectation_l,~,~,~]=MarkovChainMoments(z_grid,pi_z);
% If you look at Expectation_l you will see it is now equal to 1
Params.Expectation_l = dot(z_grid,z_prob);
%% Grids
% In the absence of idiosyncratic risk, the steady state equilibrium is given by
r_ss=1/Params.beta-1;
K_ss=((r_ss+Params.delta)/Params.alpha)^(1/(Params.alpha-1)); %The steady state capital in the absence of aggregate uncertainty.
% Set grid for asset holdings
k_grid=10*K_ss*(linspace(0,1,n_k).^3)'; % linspace ^3 puts more points near zero, where the curvature of value and policy functions is higher and where model spends more time
% Bring model into the notational conventions used by the toolkit
d_grid=0; %There is no d variable
a_grid=k_grid;
% pi_z;
% z_grid
n_d=0;
n_a=n_k;
% n_z
% Create functions to be evaluated
FnsToEvaluate.K = @(aprime,a,s) a; %We just want the aggregate assets (which is this periods state)
% Now define the functions for the General Equilibrium conditions
% Should be written as LHS of general eqm eqn minus RHS, so that the closer the value given by the function is to
% zero, the closer the general eqm condition is to holding.
GeneralEqmEqns.CapitalMarket = @(r,K,alpha,delta,Expectation_l) r-(alpha*(K^(alpha-1))*(Expectation_l^(1-alpha))-delta); %The requirement that the interest rate corresponds to the agg capital level
% Inputs can be any parameter, price, or aggregate of the FnsToEvaluate
fprintf('Grid sizes are: %i points for assets, and %i points for exogenous shock \n', n_a,n_z)
%%
DiscountFactorParamNames={'beta'};
% (a',a,s) are next-period endog state,endog state,exogenous state (Markov)
ReturnFn=@(aprime,a,s,alpha,delta,mu,r) Aiyagari1994_ReturnFn(aprime,a,s,alpha,delta,mu,r);
%%
% Use the toolkit to find the equilibrium price index
GEPriceParamNames={'r'};
% Set initial value for interest rates (Aiyagari proves that with idiosyncratic
% uncertainty, the eqm interest rate is limited above by it's steady state value
% without idiosyncratic uncertainty, that is that r<r_ss).
Params.r=0.038;
%% Test VFI
vfoptions=struct();
vfoptions.verbose = 1;
vfoptions.parallel = 0;
vfoptions.returnmatrix=1;
% Pre-compute here the return matrix
tic
if vfoptions.parallel==0 && vfoptions.returnmatrix==1
z_grid = gather(z_grid);
pi_z = gather(pi_z);
ReturnFn = zeros(n_a,n_a,n_z);
for i3=1:n_z
for i2=1:n_a
for i1=1:n_a
ReturnFn(i1,i2,i3)=Aiyagari1994_ReturnFn(a_grid(i1),a_grid(i2),z_grid(i3),Params.alpha,Params.delta,Params.mu,Params.r);
end
end
end
end
toc
% Special version for CPU
if isfield(vfoptions,'parallel') && vfoptions.parallel==0 && vfoptions.returnmatrix==0
ReturnFn=@(aprime,a,s,alpha,delta,mu,r) Aiyagari1994_ReturnFn_cpu(aprime,a,s,alpha,delta,mu,r);
end
[V,Policy]=ValueFnIter_Case1(n_d,n_a,n_z,d_grid,a_grid,z_grid,pi_z,ReturnFn,Params,DiscountFactorParamNames,[],vfoptions);
mean(V,"all")
size(Policy)
pol_aprime = squeeze(Policy(1,:,:));
pol_aprime_val = a_grid(pol_aprime);
figure
plot(a_grid,a_grid,'--')
hold on
plot(a_grid,pol_aprime_val(:,1))
hold on
plot(a_grid,pol_aprime_val(:,n_z))
