% LQRAMSM 
% 
% Compute Ramsey equilibria in a LQ economy with distorting taxation: 
% case of a Markov chain 
% The program computes allocation (consumption, leisure), tax rates, 
% revenues, and the net present value of the debt. 
% See also LQRAMSEY. 
% FRV Apr 10 1998 
%%%%%%%%%%%%%%%%% customization required %%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% 
% Gentle user: 
% The only customization required is in part (1). 
% (1) Define the parameters of the model 
% model 1: 
% 3 states 
P = [0.8 0.2 0;... % The transition matrix 
0 0.5 0.5; 
0 0 1]; 
xbar = [ 0.5 0 2.2 0 1;...% Possible states of the world: each column 
0.5 0 2.2 0 1;... % is a state of the world 
0.25 0 2.2 0 1]'; % is a state of the world 
% the rows are: [g d b s 1] 
x0=1; % index of the initial state 
Sg = [1 0 0 0 0]; % The selector matrix for the government expenditure 
Sd = [0 1 0 0 0]; % The selector matrix for the exogenous endowment 
Sb = [0 0 1 0 0]; % The selector matrix for the bliss point 
Ss = [0 0 0 1 0]; % The selector matrix for the promised coupon payments 
% Misc 
T=20; % length of simulation 
beta=1/1.05; % discount factor 
%%%%%%%%%%%%%%%%%%% This concludes the customization. Let'er rip! 
%%%%%%%%%%%%%%%% 
% (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. 
m = size(xbar,2); % Number of possible states of the world 
Sm=Sb-Sd-Ss; % a short-hand 
invP=inv(eye(size(P))-beta*P); 
a0=((1:m)==x0)*0.5*invP*((Sm*xbar)'.^2); 
b0=((1:m)==x0)*0.5*invP*(((Sb-Sd+Sg)*xbar).*((Sg+Ss)*xbar))'; 
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 
epsilon = rand(1,T); 
cumP = cumsum(P')'; % The cumulative sum for each row of P 
x(:,1) = xbar(:,x0); 
state=zeros(T,1); 
state(1)=x0; 
for t = 2:T 
state(t) = min(find(cumP(state(t-1),:) >= epsilon(t))); 
x(:,t) = xbar(:,state(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 
% compute the debt 
B=invP*(((Sb-Sc)*xbar).*((Sl-Sg)*xbar)-(Sl*xbar).^2)'; 
B=B(state,:)./p'; 
% 1-period risk-free interest rate 
R=(beta*(P(state,:)*((Sb-Sc)*xbar)')./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=p(2:T)'./(P(state(1:T-1),:)*((Sb-Sc)*xbar)'); 
Pitilde=cumsum(pi.*adjfac); 
%%%%%%% (5) plot 
set(0,'DefaultAxesColorOrder',[0 0 0],... 
'DefaultAxesLineStyleOrder','-|--|-|-.') 
figure 
subplot(2,2,1),%grid 
plot([c' g' rev']) %cons, g and rev 
v=axis; axis([1 T 0 v(4)]); 
set(gca,'FontName','Times New Roman'); 
subplot(2,2,2),%grid 
plot((1:T),g,(1:T-1),B(2:T)), %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),%grid 
plot((1:T),R-1) % interest rate 
v=axis; axis([1 T 0 v(4)]); 
set(gca,'FontName','Times New Roman'); 
% format the vertical labels for the interest rate plot 
set(gca,'YTickLabel',num2str(get(gca,'YTick')'*100,'%3.0f%%')); 
subplot(2,2,4),%grid 
plot((1:T),[rev' g'],(1:T-1),pi), 
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'); 
