% LQRAMSEY 
% 
% Compute Ramsey equilibria in a LQ economy with distorting taxation: 
% case of a stochastic difference equation 
% The program computes allocation (consumption, leisure), tax rates, 
% revenues and the net present value of the debt. 
% See also LQRAMSM. 
% FRV Apr 10 1998 
%%%%%%%%%%%%%%%%% customization required %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% Gentle user: 
% The only customization required is in part (1). 
%%%%%%% (1) Define the parameters of the model 
% Misc 
T=50; % length of simulation 
beta=1/1.05; % discount factor 
% exogenous state 
% specify a stochastic process for the state of the form 
% x(t+1)= A * x(t) + C * w(t+1) 
% where w is i.i.d. N(0,1) 
% ex1: govt spending is AR(1), state = [g ; 1] 
%mg=.35; 
%rho=.7; 
%A=[rho mg*(1-rho); 0 1]; 
%C=[mg/10*sqrt(1-rho^2);0]; 
%% These selector matrices are all 1*nx, where nx is the size of the state: 
%Sg = [1 0]; % The selector matrix for the government expenditure 
%Sd = [0 0]; % The selector matrix for the exogenous endowment 
%Sb = [0 2.135]; % The selector matrix for the bliss point 
%Ss = [0 0]; % The selector matrix for the promised coupon payments 
% ex2: seasonal 
rho=.95; 
mg=.35; 
A=[0 0 0 rho mg*(1-rho); 
1 0 0 0 0 ; 
0 1 0 0 0 ; 
0 0 1 0 0; 
0 0 0 0 1]; 
C=[mg/8*sqrt(1-rho^2);0;0;0;0]; 
Sg = [1 0 0 0 0]; 
Sd = [0 0 0 0 0]; 
Sb = [0 0 0 0 2.135]; % this is chosen so that (Sc+Sg)*x0=1 
Ss = [0 0 0 0 0]; 
%%%%%%%%%%%%%%%%%%% This concludes the customization. Let'er rip! 
%%%%%%%%%%%%%%%% 
%%%%%%% (1) generate initial condition 
[nx,nx]=size(A); 
if ~exist('x0') 
x0=null(eye(nx)-A); 
if (x0(nx)<0) 
x0=-x0; 
end 
x0=x0./x0(nx); 
end 
%%%%%%% (2) solve for the Lagrange multiplier on the government BC 
% we actually solve for mu=lambda/(1+2*lambda) 
% mu is the solution to a quadratic equation a(mu^2-mu)+b=0 
% where a and b are expected discounted sums of quadratic forms 
% of the state; we use dlyap to compute those sums quickly. 
Sm=Sb-Sd-Ss; % a short-hand 
Qa=dlyap(sqrt(beta)*A',-0.5*Sm'*Sm); 
qa=trace(C'*Qa*C)*beta/(1-beta); 
Qb=dlyap(sqrt(beta)*A',-0.5*(Sb-Sd+Sg)'*(Sg+Ss)); 
qb=trace(C'*Qb*C)*beta/(1-beta); 
a0=x0'*Qa*x0+qa; 
b0=x0'*Qb*x0+qb; 
if a0<0 % normalize 
a0=-a0; b0=-b0; 
end 
disc=a0^2-4*a0*b0; 
% there may be no solution: 
if ( disc < 0 ) 
disp('There is no Ramsey equilibrium.'); 
disp('(Hint: you probably set government spending too high;'); 
disp(' elect a Republican Congress and start over.)'); 
return 
end 
mu=0.5*(a0-sqrt(disc))/a0; 
% or it may be of the wrong sign: 
if ( mu*(0.5-mu) < 0 ) 
disp('negative multiplier on the government budget constraint!'); 
disp('(Hint: you probably set government spending too low;'); 
disp(' elect a Democratic Congress and start over.)'); 
return 
end 
%%%%%%% (3a) solve for the allocation 
Sc=0.5*(Sb+Sd-Sg-mu*Sm); 
Sl=0.5*(Sb-Sd+Sg-mu*Sm); 
Stau1=2*mu*Sm; 
Stau2=Sb-Sd+Sg+mu*Sm; 
%%%%%%% (4) run the simulation 
[nx,nw]=size(C); 
x=zeros(nx,T); 
w=randn(nw,T); % draw the shocks 
x(:,1)=x0; % initialize the state 
for t=2:T % compute x recursively 
x(:,t)=(A*x(:,t-1) + C*w(:,t)); 
end 
% exogenous stuff: 
g = Sg*x; % govmint spending 
d = Sd*x; % endowment 
b = Sb*x; % bliss point 
s = Ss*x; % coupon payment on existing debt 
% endogenous stuff: 
c=Sc*x; % consumption 
l=Sl*x; % labor 
p=(Sb-Sc)*x; % price 
tau=(Stau1*x)./(Stau2*x); % taxes 
rev=tau.*l; % revenues 
QB=dlyap(sqrt(beta)*A',-((Sb-Sc)'*(Sl-Sg)-Sl'*Sl)); 
qB=trace(C'*QB*C)*beta/(1-beta); 
B=(x'*QB*x+qB)/p; % debt 
% 1-period risk-free interest rate 
R=((beta*(Sb-Sc)*A*x)./p)'.^(-1); 
pi=B(2:T)-R(1:T-1).*(B(1:T-1)-rev(1:T-1)'+g(1:T-1)'); 
% risk-adjusted martingale 
adjfac=(((Sb-Sc)*x(:,2:T))./((Sb-Sc)*A*x(:,1:T-1)))'; 
Pitilde=cumsum(pi.*adjfac); 
%%%%%%% (5) plot 
set(0,'DefaultAxesColorOrder',[0 0 0],... 
'DefaultAxesLineStyleOrder','-|--|-|-.') 
figure 
subplot(2,2,1), 
plot([c' g' rev']),grid %cons, g and rev 
v=axis; axis([1 T 0 v(4)]); 
set(gca,'FontName','Times New Roman'); 
subplot(2,2,2) 
plot((1:T),g,(1:T-1),B(2:T)),grid %g(t) and B(t+1) 
v=axis; axis([1 T v(3) v(4)]); 
set(gca,'FontName','Times New Roman'); 
subplot(2,2,3) 
plot((1:T),R-1),grid % interest rate 
v=axis; axis([1 T v(3) v(4)]); 
set(gca,'FontName','Times New Roman'); 
% format the vertical labels for the interest rate plot 
yticks=get(gca,'YTick')'; 
set(gca,'YTickLabel',num2str(yticks*100,'%4.1f%%')); 
subplot(2,2,4) 
plot((1:T),[rev' g'],(1:T-1),pi), grid 
v=axis; axis([1 T v(3) v(4)]); 
set(gca,'FontName','Times New Roman'); 
figure 
subplot(2,2,1), 
plot((2:T),adjfac),grid 
v=axis; axis([1 T v(3) v(4)]); 
set(gca,'FontName','Times New Roman'); 
subplot(2,2,2),grid 
plot((2:T),Pitilde) 
v=axis; axis([1 T v(3) v(4)]); 
set(gca,'FontName','Times New Roman'); 
% th-th-th-that's all, folks! 
return 
% labels to place on the graphs 
gtext('\pi\it(t+1)'); 
gtext('\it g(t)'); 
gtext('\it g(t)'); 
gtext('\it g(t)'); 
gtext('\it c(t)'); 
gtext('\it\tau(t)l(t)'); 
gtext('\it\tau(t)l(t)'); 
h=gtext('\it B(t+1)'); 
h=gtext('\it R(t)-1'); 
