Please enable javascripts for pretty navigation. Here are ugly links to make it usable:

Example file of how to use adaptive sparse grid for HJB Equation

by SeHyoun Ahn, July 2016

Contents

% Requires: <generate_spgrid.m>
%
%           <gen_diff_op.m>
%
%           <grid_evaluation.m>
%

Set up Grids

n_dim = 2;                 % Dimensionality of the problem
n_levels = 4;              % Level of the starting sparse grid
I = 100;                   % Number of asset grid points for full grid
J = 40;                    % Number of skill grid points for full grid

% Generate sparse grid
[spgrid,spspace] = generate_spgrid(n_dim,n_levels);


% Setup Change of basis from Hierachical to Nodal
hier_to_nodal = grid_evaluation(spgrid,spgrid,spspace);

% Generate Difference Operators
[forward_diff_x,backward_diff_x] = gen_diff_op(1,spgrid,spspace,min(spspace(:,1)));
[forward_diff_y,backward_diff_y] = gen_diff_op(2,spgrid,spspace,min(spspace(:,2)));

% Second Derivative Operator
% More tests are needed to check which stencil is the best for second order
% Operators
second_deriv_y = forward_diff_y*(hier_to_nodal\backward_diff_y) ...
  + backward_diff_y*(hier_to_nodal\forward_diff_y);
second_deriv_y(spgrid(:,2)==1,:) = -backward_diff_y(spgrid(:,2)==1,:);
second_deriv_y(spgrid(:,2)==0,:) = forward_diff_y(spgrid(:,2)==0,:);

% Integration Weights
% This does not get used here
int_weight = (prod(spspace,2)./2.^sum(spspace==1,2))';

Set up Model

Parameters

ga = 2;             % CRRA utility with parameter gamma
alpha = 0.35;       % Production function F = K^alpha * L^(1-alpha)
delta = 0.1;        % Capital depreciation
zmean = 1.0;        % mean O-U process (in levels). This parameter has to be adjusted to ensure that the mean of z (truncated gaussian) is 1.
sig2 = (0.10)^2;    % sigma^2 O-U
Corr = exp(-0.3);   % persistence -log(Corr)  O-U
rho = 0.05;         % discount rate
TFP = 1;

K = 3.8;            % initial aggregate capital. It is important to guess a value close to the solution for the algorithm to converge

zmin = 0.5;   % Range z
zmax = 1.5;
amin = -1;    % borrowing constraint
amax = 30;    % range a

%simulation parameters
maxit  = 100;       %maximum number of iterations in the HJB loop
maxitK = 100;       %maximum number of iterations in the K loop
crit = 10^(-10);    %criterion HJB loop
critK = 1e-7;       %criterion K loop
Delta = 1000;       %delta in HJB algorithm

% Ornstein-Uhlenbeck
the = -log(Corr);
Var = sig2/(2*the);

% Evalute asset and skill values for grid points
aa = spgrid(:,1)*(amax-amin)+amin;
zz = spgrid(:,2)*(zmax-zmin)+zmin;
mu = the*(zmean-zz);
s2 = sig2;

% Differential operator for skill
Aswitch = bsxfun(@times,mu.*(mu>=0),forward_diff_y)/(zmax-zmin) ...
  + bsxfun(@times,mu.*(mu<0),backward_diff_y)/(zmax-zmin) ...
  + s2/2*second_deriv_y/(zmax-zmin)^2;

% Initial guess
r = alpha     * TFP * K^(alpha-1) -delta; %interest rates
w = (1-alpha) * TFP * K^(alpha);          %wages
v0 = (w.*zz + r.*aa).^(1-ga)/(1-ga)/rho;
v = hier_to_nodal\v0;
dist = zeros(1,maxit);
V=v;

number_of_points=0;

Iterate to Solve HJB Equation Adding Nodes to Sparse Grids as Needed

Note that the Iteration Part of the Code is the same except for difference operator

while true
  for n=1:maxit
    Vaf = forward_diff_x*V/(amax-amin);
    Vab = backward_diff_x*V/(amax-amin);

    c0 = w*zz + r.*aa;
    cf = Vaf.^(-1/ga);
    sf = w*zz + r.*aa - cf;
    cb = Vab.^(-1/ga);
    sb = w*zz + r.*aa - cb;

    If = (sf >  0).*(spgrid(:,1)<1);
    Ib = (sb < 0).*(If==0).*(spgrid(:,1)>0);
    I0 = (1-If-Ib);

    c = cf.*If + cb.*Ib + c0.*I0;
    u = c.^(1-ga)/(1-ga);

    AA = bsxfun(@times,If.*sf,forward_diff_x)/(amax-amin) + bsxfun(@times,Ib.*sb,backward_diff_x)/(amax-amin);

    B = (1/Delta + rho)*hier_to_nodal - AA - Aswitch;

    b = u + hier_to_nodal*V/Delta;

    V_new = B\b;

    Vchange = V_new - V;
    V = V_new;

    dist = max(max(abs(Vchange)));
    if dist<crit
      disp('Value Function Converged, Iteration = ')
      disp(n)
      break
    end
  end

  % Plot Value function given the grid
  figure(1);
  scatter3(aa,zz,hier_to_nodal*V);
  drawnow;

  V_nodal = hier_to_nodal*V;
Value Function Converged, Iteration = 
     9

Add new nodes to sparse grid if needed

  % Find Position to Expand
  loc = abs(V)>1e-4*(max(V_nodal)-min(V_nodal));

  % Compute New Grid Points from active grid points
  [new_points,new_spaces] = gen_new_points(spgrid(loc,:),spspace(loc,:));

  % Find Positions to Compress Away
  loc_keep = abs(V)>8e-5*(max(V_nodal)-min(V_nodal));
  [xx,loc] = unique([spgrid(loc_keep,:);new_points],'rows');

  % Break if no new grid points are added
  if number_of_points >= size(xx,1)
    break;
  end
  number_of_points = size(xx,1);

  % Update interpolation operator
  x_eval = grid_evaluation(xx,spgrid,spspace);
  v_nodal = x_eval*V;
  spgrid = xx;
  spspace = [spspace(loc_keep,:);new_spaces];
  spspace = spspace(loc,:);
  hier_to_nodal = grid_evaluation(xx,spgrid,spspace);

  % Correct the basis function for the new bases
  V = hier_to_nodal\v_nodal;

  % Update Difference Operators for the new grid
  [forward_diff_x,backward_diff_x] = gen_diff_op(1,spgrid,spspace,min(spspace(:,1)));
  [forward_diff_y,backward_diff_y] = gen_diff_op(2,spgrid,spspace,min(spspace(:,2)));
  second_deriv_y = forward_diff_y*(hier_to_nodal\backward_diff_y) ...
    + backward_diff_y*(hier_to_nodal\forward_diff_y);
  second_deriv_y(spgrid(:,2)==1,:) = -backward_diff_y(spgrid(:,2)==1,:);
  second_deriv_y(spgrid(:,2)==0,:) = forward_diff_y(spgrid(:,2)==0,:);

  % Update Nodal Values for new grid
  aa = spgrid(:,1)*(amax-amin)+amin;
  zz = spgrid(:,2)*(zmax-zmin)+zmin;
  mu = the*(zmean-zz);
  s2 = sig2;

  Aswitch = bsxfun(@times,mu.*(mu>=0),forward_diff_y)/(zmax-zmin) ...
    + bsxfun(@times,mu.*(mu<0),backward_diff_y)/(zmax-zmin) ...
    + s2/2*second_deriv_y/(zmax-zmin)^2;
end
  
Value Function Converged, Iteration = 
     7

Value Function Converged, Iteration = 
     8

Value Function Converged, Iteration = 
     8

Value Function Converged, Iteration = 
     9

Value Function Converged, Iteration = 
     8

Value Function Converged, Iteration = 
     8

Value Function Converged, Iteration = 
     7

Value Function Converged, Iteration = 
     6

Value Function Converged, Iteration = 
     7

figure(2);
scatter3(aa,zz,hier_to_nodal*V);
V_sparse = V;

Comparison Plot to Full Grid

The following code is adapted from Ben Moll's codes available at https://www.princeton.edu/~moll/HACTproject.htm

a = linspace(amin,amax,I)';    % wealth vector
da = (amax-amin)/(I-1);
z = linspace(zmin,zmax,J);     % productivity vector
dz = (zmax-zmin)/(J-1);
dz2 = dz^2;
aa = a*ones(1,J);
zz = ones(I,1)*z;

mu = the*(zmean - z);        % Drift
s2 = sig2.*ones(1,J);        % Variance

% Preallocation
Vaf = zeros(I,J);
Vab = zeros(I,J);
Vzf = zeros(I,J);
Vzb = zeros(I,J);
Vzz = zeros(I,J);
c = zeros(I,J);

yy = - s2/dz2 - mu/dz;
chi =  s2/(2*dz2);
zeta = mu/dz + s2/(2*dz2);

updiag=zeros(I*J,1); %This is necessary because of the peculiar way spdiags is defined.
updiag(I+1:I*J)=reshape(repmat(zeta(1:J-1),I,1),I*(J-1),1);

centdiag=zeros(I*J,1);
centdiag(1:I)=repmat(chi(1)+yy(1),I,1);
centdiag(I+1:I*(J-1))=reshape(repmat(yy(2:J-1),I,1),I*(J-2),1);
centdiag(I*J-I+1:I*J)=repmat(yy(J)+zeta(J),I,1);

lowdiag=reshape(repmat(chi(2:J),I,1),I*(J-1),1);

Aswitch=spdiags(centdiag,0,I*J,I*J)+spdiags(lowdiag,-I,I*J,I*J)+spdiags(updiag,I,I*J,I*J);

% Initial Guess
r = alpha     * TFP * K^(alpha-1) -delta; %interest rates
w = (1-alpha) * TFP * K^(alpha);          %wages
v0 = (w*zz + r.*aa).^(1-ga)/(1-ga)/rho;
v = v0;

for n=1:maxit
  V = v;

  % forward difference
  Vaf(1:I-1,:) = (V(2:I,:)-V(1:I-1,:))/da;
  Vaf(I,:) = (w*z + r.*amax).^(-ga); %will never be used, but impose state constraint a<=amax just in case

  % backward difference
  Vab(2:I,:) = (V(2:I,:)-V(1:I-1,:))/da;
  Vab(1,:) = (w*z + r.*amin).^(-ga);  %state constraint boundary condition

  %consumption and savings with forward difference
  cf = Vaf.^(-1/ga);
  sf = w*zz + r.*aa - cf;
  %consumption and savings with backward difference
  cb = Vab.^(-1/ga);
  sb = w*zz + r.*aa - cb;
  %consumption and derivative of value function at steady state
  c0 = w*zz + r.*aa;
  Va0 = c0.^(-ga);

  If = sf > 0; %positive drift --> forward difference
  Ib = sb < 0; %negative drift --> backward difference
  I0 = (1-If-Ib); %at steady state

  Va_Upwind = Vaf.*If + Vab.*Ib + Va0.*I0; %important to include third term

  c = Va_Upwind.^(-1/ga);
  u = c.^(1-ga)/(1-ga);

  X = - min(sb,0)/da;
  Y = - max(sf,0)/da + min(sb,0)/da;
  Z = max(sf,0)/da;

  Z(I,:)=0;
  X(1,:)=0;
  updiag=[0;reshape(Z,I*J,1)];
  centdiag=reshape(Y,I*J,1);
  lowdiag=reshape(X,I*J,1);

  AA=spdiags(centdiag,0,I*J,I*J)+spdiags([updiag;0],1,I*J,I*J)+spdiags(lowdiag(2:I*J),-1,I*J,I*J);
  A = AA + Aswitch;
  B = (1/Delta + rho)*speye(I*J) - A;

  u_stacked = reshape(u,I*J,1);
  V_stacked = reshape(V,I*J,1);

  b = u_stacked + V_stacked/Delta;

  V_stacked = B\b;

  V = reshape(V_stacked,I,J);

  Vchange = V - v;
  v = V;

  dist = max(max(abs(Vchange)));
  if dist<crit
    disp('Value Function Converged with Full Grid, Iteration = ');
    disp(n);
    break;
  end
end

h = figure(2);
hold on;
surf(a,z,V');
title('Comparison to Uniform Grid');
gif_maker_rotate('comparison.gif',h);
Value Function Converged with Full Grid, Iteration = 
     8