Coleman (1991): A production economy with an income tax

Coleman (1991) studies a variation of the RBC model with distortionary taxes. The author proves the existence and uniqueness of recursive equilibrium and charaterizes its properties in this economy. In addition, the author provides an algorithm to compute recursive equilibrium using policy function iteration and proves its convergence. This is one of the rare examples in which the convergence of the policy function iteration can be proved analytically.

Students in the advanced macroeconomics class (Econ 608) at Georgetown University in Spring 2020 (Tae-Hun Chang, Chen Jin, Menghan Li, Yingqi Xu) implented this algorithm in the GDSGE gmod file below.

The gmod File

Coleman1991.gmod

  1
  2
  3
  4
  5
  6
  7
  8
  9
 10
 11
 12
 13
 14
 15
 16
 17
 18
 19
 20
 21
 22
 23
 24
 25
 26
 27
 28
 29
 30
 31
 32
 33
 34
 35
 36
 37
 38
 39
 40
 41
 42
 43
 44
 45
 46
 47
 48
 49
 50
 51
 52
 53
 54
 55
 56
 57
 58
 59
 60
 61
 62
 63
 64
 65
 66
 67
 68
 69
 70
 71
 72
 73
 74
 75
 76
 77
 78
 79
 80
 81
 82
 83
 84
 85
 86
 87
 88
 89
 90
 91
 92
 93
 94
 95
 96
 97
 98
 99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
% Parameters
parameters beta sigma alpha delta K_min K_max rho Imin tau;

beta = 0.99; % dicount factor
sigma = 2; % risk aversion
alpha = 0.36; % capital share
delta = 0.025; % depreciation rate
rho = 0.90; %persistent of productivity shocks
phi = 0.975; %investment irreversibility constraint
tau = 0.2; %linear tax income

TolEq = 1e-6;
SaveFreq = 50;
SimuSaveFreq = 1000;
SimuPrintFreq = 1000;

NumThreads = feature('NumCores');

var_shock e;

% Shocks
shock_num = 2;
shock_trans = [1/2.0,1/2.0;1/2.0,1/2.0];

e = [-0.01,0.01]*sqrt(1-rho^2);

% Steady state
Kss = (alpha/(1/beta - 1 + delta))^(1/(1-alpha));
Iss = delta*Kss;
Imin = phi*Iss;

K_min = 0.5*Kss;
K_max = 1.5*Kss;    

% States
var_state z K;

z_pts = 21;
z = linspace(0.9,1.1,z_pts);

K_pts = 101;
K = exp(linspace(log(K_min),log(K_max),K_pts));

var_interp Euler_interp muc_interp;

% Last period
var_policy_init c;

inbound_init c 1e-6 100;

var_aux_init Euler mu muc;

model_init;
  
  Y = z*(K^alpha);
  
  resid = 1 - c/(Y-Imin);
  
  mu = 1;

  muc = mu*(c^-sigma);

  Euler = (c^-sigma)*(1-delta+(1-tau)*z*alpha*K^(alpha-1));
  
  equations;
     resid;		
  end;
end;

initial Euler_interp Euler;

Euler_interp = Euler;

initial muc_interp muc;

muc_interp = muc;

% Endogenous variables, bounds, and initial values

var_policy c invst mu;

inbound c 1e-6 100 adaptive(2);
inbound invst 0 50 adaptive(1.5);
inbound mu 0 2 adaptive(1.5);

% Other equilbirium variables

var_aux Y Euler K_next muc z_next[2] Inv;

model;

  Inv = invst + Imin;
     
  K_next = Inv + (1-delta)*K;

  Y = z*(K^alpha);
  
  resid  = 1 - (Inv+c)/Y;

  % transition to z_next

  z_next' = exp(rho*log(z)+e');
   
  Euler_future' = Euler_interp'(z_next',K_next);
  
  EEuler_future = GNDSGE_EXPECT{Euler_future'};
   
  muc_future' = muc_interp'(z_next',K_next);

  Emuc_future = GDSGE_EXPECT{muc_future'};
   
  Euler = (c^-sigma)*(1-delta+(1-tau)*z*alpha*K^(alpha-1));

  % First order conditions for firms
  
  foc = 1-mu - (beta*EEuler_future-beta*(1-delta)*Emuc_future)/(c^-sigma);
  
  comp_slack = mu*invst;

  muc = mu*(c^-sigma);

  equations;
      foc;
      resid;
      comp_slack; 
  end;

end;

simulate;

  num_periods = 15000;
  num_samples = 100;
  
  initial K Kss;
  initial z 1;

  initial shock 2;
  
  var_simu z K Y c Inv;

  K' = K_next;
  z' = z_next';

end;