Is there a MATLAB equivalent to QuantEcon.LinInterp?
I am rewriting the value function iteration code from the Julia lessons on single agent models into MATLAB. I am hoping to better understand how VFI works and then use it to solve an asset pricing problem slightly more complicated than the income fluctuation problem.
The MATLAB equivalent would probably be interp1q (https://www.mathworks.com/help/matlab/ref/interp1q.html).
You might find some of our workshop material helpful for understanding VFI and using it to solve models. See this repository and, in particular, this very nice notebook by Victoria Gregory.
I discovered the answer independently of Lutz Hendricks, who also notes that griddedInterpolant enables one to write MATLAB code with a solution similar to the Julia version on QuantEcon.
His code is much prettier than mine and can be found here: https://github.com/hendri54/Econ821/blob/master/matlab_examples/growth_model_821.m
Below is the code I ended up writing:
Bellman Operator Function:
function [sigma, Tw] = stochgrowth1_bellman_operator(grid_y,... % Considered y-values (discretized).
w_0,... % Initial guess for associated value function.
alpha,... % Elasticity of capital.
betta,... % Discount rate.
mu,... % Mean of shock process.
S,... % Standard error of shock process.
shsize,... % Number of shock draws in Monte Carlo integral (optional).
tol,... % Optimizer tolerance (optional).
u,... % Utility function (optional).
f) % Production function (optional)
% Bellman operator function for stochastic growth model with one state
% variable.
% Check number of inputs
if nargin > 10
error('stochgrowth_bellman_operator:TooManyInputs', ...
'requires at most 4 optional inputs');
end
% Fill in unset optional values.
switch nargin
case 6
shsize = 250;
tol = 1e-6;
u = @(c) log(c);
f = @(k) k^alpha;
case 7
tol = 1e-6;
u = @(c) log(c);
f = @(k) k^alpha;
case 8
u = @(c) log(c);
f = @(k) k^alpha;
case 9
f = @(k) k^alpha;
end
% Set tolerance for optimizer.
optS = optimset('fminbnd');
optS.TolX = tol;
% Output of the Bellman operator for each k grid point.
sigma = NaN(size(w_0));
Tw = sigma;
% Shock process
shocks = exp(mu + S * randn(shsize,1));
% Interpolation object.
w_func = griddedInterpolant(grid_y, w_0, 'linear');
for i_y = 1:length(grid_y)
% Given level of income.
y = grid_y(i_y);
% Find acceptable range for c.
cMin = 1e-6;
cMax = y - cMin;
% Solve for extremal value of functional,
% i.e. V(y_t) = max( u(c) + betta*E[V(y_t+1)] ),
% given a particular level of income.
objective = @(c) -u(c) - betta * mean(w_func(f(y-c) * shocks));
[X, FVAL] = fminbnd(objective, cMin, cMax);
% Return values.
sigma(i_y) = X; % Policy function.
Tw(i_y) = -FVAL; % Value function.
end % End for-loop.
end % End bellman_operator function.
Main Script:
%% stochastic growth model
% In theory, the algorithm is as follows:
% 1. Begin with a function w() -- an initial condition
% 2. Solving the Bellman operator, obtain the function Tw()
% 3. Unless some stopping condition is satisfied, set w = Tw and return
% to step 2.
%% housekeeping
close all
clear all
%% parameters
alpha = 0.40; % Elasticity of capital.
betta = 0.96; % Discount rate.
mu = 0.00; % Mean of shock process.
S = 0.10; % Standard error of shock process.
niter = 35; % Number of iterations to find (converged) value function.
%% discretize space
grid_max = 4; % largest grid point
grid_size = 200; % number of grid points
grid_y = linspace(1e-5, grid_max, grid_size); % income level state-space under consideration
%% Test stochgrowth_bellman_operator() by replicating analytic solution numerically.
% Initial guess for V -- analytic solution.
c1 = log(1 - alpha*betta) / (1 - betta);
c2 = (mu + alpha * log(alpha*betta)) / (1 - alpha);
c3 = 1 / (1 - betta);
c4 = 1 / (1 - alpha * betta);
w_0 = v_star(grid_y,c1,c2,c3,c4);
% Find extremal of functional given state-space for income and other
% parameter values.
[sigma, Tw] = stochgrowth1_bellman_operator(grid_y,w_0,alpha,betta,mu,S);
% Plot results.
figure()
ax1 = subplot(2,2,1);
plot(ax1,grid_y,sigma)
ylim([0 2.5])
title('Policy function')
xlabel('Income')
ylabel('Consumption')
grid on
ax3 = subplot(2,2,3);
plot(ax3, grid_y,Tw)
ylim([-40 10])
title('Value function')
xlabel('Income')
ylabel('Value')
grid on
%% solve for optimal value function and policy function, given naive starting value.
% Discretize state-space for income.
grid_max = 4; % largest grid point
grid_size = 200; % number of grid points
grid_y = linspace(1e-5, grid_max, grid_size); % income level state-space under consideration
% guess initial condition for V().
w_0 = 5 * log(grid_y);
% Create blank storage for policy function and value function
Sigma = NaN(grid_size,niter);
V = Sigma;
% Iterate Bellman operator, updating guess for value function
for iter = 1:niter
% Solve Bellman operator.
[X, FVAL] = stochgrowth1_bellman_operator(grid_y,w_0,alpha,betta,mu,S);
% Get values.
Sigma(:,iter) = X; % Policy function.
V(:,iter) = FVAL; % Value function.
% Update initial condition.
w_0 = FVAL;
end
% Graph properties.
% COLORS = linspecer(niter); % see: https://www.mathworks.com/matlabcentral/fileexchange/42673-beautiful-and-distinguishable-line-colors-+-colormap
% Plot graphs.
ax2 = subplot(2,2,2);
% set(ax2, 'ColorOrder', COLORS)
hold on
plot(ax2, grid_y,Sigma)
ylim([0 2.5])
title('Policy function')
xlabel('Income')
ylabel('Consumption')
grid on
ax4 = subplot(2,2,4);
% set(ax4, 'ColorOrder', COLORS)
hold on
plot(grid_y,V)
ylim([-40 10])
title('Value function')
xlabel('Income')
ylabel('Value')
grid on