%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
% solves dynamic ramsey taxation problem         %
%                                                %
% florian scheuer, 2007                          %
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

clear all

global sigma theta alpha delta beta A

tic

%% parameters

sigma=.5;
beta=.947;
theta=.36;
delta=.08;
alpha=.4;
A=1;

% taxes
tau=0;
kappa=0;
kappagrid=200;
kappa_vec=[-.5:1/kappagrid:.5]';
gridn=length(kappa_vec);
rebate=1; %1 if lump sum rebate of revenue, 0 if revenue spent on govt consumption

%% compute steady state as function of taxes

R=1/beta;
r_vec=(R-1)./(1-kappa_vec)+delta;
kl_vec=(alpha*A./r_vec).^(1/(1-alpha));
w_vec=(1-alpha).*A.*kl_vec.^alpha;

if rebate==1
    cl_vec=A.*kl_vec.^alpha-delta.*kl_vec;
else
    gl_vec=tau.*w_vec+kappa.*(r_vec-delta).*kl_vec;    
    cl_vec=A.*kl_vec.^alpha-delta.*kl_vec-gl_vec;
end

l_vec=(theta.*(1-tau).*w_vec)./((1-theta).*cl_vec+theta.*(1-tau).*w_vec);
c_vec=cl_vec.*l_vec;
k_vec=kl_vec.*l_vec;
y_vec=A.*k_vec.^alpha.*l_vec.^(1-alpha);
u_vec=u(c_vec,l_vec);   
v_vec=u_vec/(1-beta);

% lambda
if sigma==1
    lambda_vec=(exp((1-beta)*v_vec(kappagrid/2+1))./(1-l_vec).^(1-theta)).^(1/theta)./c_vec-1;
else
    lambda_vec=(((1-beta)*(1-sigma)*v_vec(kappagrid/2+1)).^(1/(1-sigma))./(1-l_vec).^(1-theta)).^(1/theta)./c_vec-1;
end

% plot

figure(1); subplot(3,1,1); plot(kappa_vec',[c_vec'; l_vec'; y_vec']); title('steady state c, l, y as function of kappa'); grid on
figure(1); subplot(3,1,2); plot(kappa_vec',k_vec'); title('steady state k as function of kappa'); grid on
figure(1); subplot(3,1,3); plot(kappa_vec',lambda_vec'); title('lambda as function of kappa'); grid on

%% value function iteration without taxes

% set capital state space interval as fractions of steady state
klow    = min(k_vec); khigh = max(k_vec);

% grid size and vector of ones
% gridn   = kappagrid;
e       = ones(gridn,1);

% create grid for capital
%k = (klow: (khigh-klow)/(gridn-2) : khigh)';
kss = k_vec(kappagrid/2+1);
css=c_vec(kappagrid/2+1);
k = k_vec;      

k_mat=e*k';
kprime_mat=k*e';

% minimal l so that the combination k, k' is feasible
lmin = min(((max(0,kprime_mat-(1-delta).*k_mat))./(A.*k_mat.^alpha)).^(1/(1-alpha)),1);

% find the optimal l. create matrix U(k,k')
U=zeros(gridn,gridn);
for i=1:gridn
    for j=1:gridn
        if lmin(i,j)==1;
            U(i,j)=-1e20;
        else
            lopt=fminbnd(@(l) -u(pf(k_mat(i,j),l)-kprime_mat(i,j),l),lmin(i,j),1);
            U(i,j)=u(pf(k_mat(i,j),lopt)-kprime_mat(i,j),lopt);
        end
    end
end

% initialize iteration
crit    = 1; iter = 1; 
l_init=.5;
v  = u(pf(k,l_init) - k,l_init)/(1-beta);  % initial guess (value of keeping k and c constant)
%v =0*v;                    % alternative initial guess

% iterate on Bellman equation
while crit > 1e-10
    
    vnew = max( U + beta*v(:,iter)*e' )';

    % compute normalized criterion (scale is in consumption per period)
    critnew = (1-beta)*max(abs(vnew - v(:,iter)))/css^(-sigma+1);
        
    % display iteration information
    % note: contraction mapping implies that critnew/crit < beta
    [iter , critnew, critnew/crit ]
    crit = critnew; iter = iter +1; v(: , iter)  = vnew;
    
end

% compute policy functions
[v, policykindex]= max(U + beta*v(:,iter)*e');
policyk = k(policykindex);
policyl=zeros(length(k),1);
policyc=zeros(length(k),1);
for i=1:length(k)
    policyl(i,1)=fminbnd(@(l) -u(pf(k(i,1),l)-policyk(i,1),l),lmin(policykindex(i),i),1);
    policyc(i,1)=pf(k(i,1),policyl(i,1))-policyk(i,1);
end
v=v';

% lambda
if sigma==1
    lambda_vec2=(exp((1-beta)*v)./(1-l_vec).^(1-theta)).^(1/theta)./c_vec-1;
else
    lambda_vec2=(((1-beta)*(1-sigma)*v).^(1/(1-sigma))./(1-l_vec).^(1-theta)).^(1/theta)./c_vec-1;
end


figure(2); subplot(2,1,1); plot(k',v'); title('value without taxes as function of k(kappa)'); grid on
figure(2); subplot(2,1,2); plot(kappa_vec',lambda_vec2'); title('lambda as function of kappa'); grid on


%% compute ramsey optimum

mu=0.2092;
crit_mu=1;
crit_iter=1;
T=500;

kappa=.35;
tau=.35;

% initial steady state values
R0=1/beta;
r0=(R0-1)./(1-kappa)+delta;
kl0=(alpha*A./r0).^(1/(1-alpha));
w0=(1-alpha).*A.*kl0.^alpha;
gl0=tau.*w0+kappa.*(r0-delta).*kl0;    
cl0=A.*kl0.^alpha-delta.*kl0-gl0;

l0=(theta.*(1-tau).*w0)./((1-theta).*cl0+theta.*(1-tau).*w0);
c0=cl0.*l0;
k0=kl0.*l0;
g0=gl0.*l0;

% time 1 problem

% minimal l so that the combination k, k' is feasible
lmin = min(((max(0,g0+kprime_mat-(1-delta).*k_mat))./(A.*k_mat.^alpha)).^(1/(1-alpha)),1);


% iterate on mu
while crit_mu>.01 & crit_iter<10

% find the optimal l. create matrix W1(k,k')
W1=zeros(gridn,gridn);
for i=1:gridn
    for j=1:gridn
        if lmin(i,j)==1;
            W1(i,j)=-1e50;
        else
            lopt=fminbnd(@(l) -w1(pf(k_mat(i,j),l)-kprime_mat(i,j)-g0,l,mu),lmin(i,j),1);
            W1(i,j)=w1(pf(k_mat(i,j),lopt)-kprime_mat(i,j)-g0,lopt,mu);
        end
    end
end

% initialize iteration
crit    = 1; iter = 1; 
l_init=.5;
v1  = w1(pf(k,l_init) - k-g0,l_init,mu)/(1-beta);  % initial guess (value of keeping k and c constant)
%v =0*v;                    % alternative initial guess

% iterate on Bellman equation
while crit > 1e-10
    
    v1new = max( W1 + beta*v1(:,iter)*e' )';

    % compute normalized criterion (scale is in consumption per period)
    critnew = (1-beta)*max(abs(v1new - v1(:,iter)))/c0^(-sigma+1);
        
    % display iteration information
    % note: contraction mapping implies that critnew/crit < beta
    [iter , critnew, critnew/crit ]
    crit = critnew; iter = iter +1; v1(: , iter)  = v1new;
    
end

% compute policy functions
[v1, policykindex1]= max(W1 + beta*v1(:,iter)*e');
policyk1 = k(policykindex1);
policyl1=zeros(length(k),1);
policyc1=zeros(length(k),1);
for i=1:length(k)
    policyl1(i,1)=fminbnd(@(l) -w1(pf(k(i,1),l)-policyk1(i,1)-g0,l,mu),lmin(policykindex1(i),i),1);
    policyc1(i,1)=pf(k(i,1),policyl1(i,1))-policyk1(i,1)-g0;
end
v1=v1';

% time 0 problem

% find the optimal l. create matrix W0(k,k')
W0=zeros(gridn,gridn);
for i=1:gridn
    for j=1:gridn
        if lmin(i,j)==1;
            W0(i,j)=-1e20;
        else
            lopt=fminbnd(@(l) -w00(pf(k_mat(i,j),l)-kprime_mat(i,j)-g0,l,mu,k_mat(i,j)),lmin(i,j),1);
            W0(i,j)=w00(pf(k_mat(i,j),lopt)-kprime_mat(i,j)-g0,lopt,mu,k_mat(i,j));
        end
    end
end

% compute v0 and policy functions
[v0,policykindex0]=max(W0 + beta*v1*e');
v0=v0';
policyk0=k(policykindex0);
policyl0=zeros(length(k),1);
policyc0=zeros(length(k),1);
for i=1:length(k)
    policyl0(i,1)=fminbnd(@(l) -w00(pf(k(i,1),l)-policyk0(i,1)-g0,l,mu,k(i,1)),lmin(policykindex0(i),i),1);
    policyc0(i,1)=pf(k(i,1),policyl0(i,1))-policyk0(i,1)-g0;
end

% compute time series for c,l,k

k_t=zeros(T+1,1);
l_t=zeros(T,1);
c_t=zeros(T,1);

k_t(1,1)=k0;
index_k0=dsearchn(k,k0);
c_t(1,1)=policyc0(index_k0);
l_t(1,1)=policyl0(index_k0);
k_t(2,1)=policyk0(index_k0);

for t=2:T
    index_k=dsearchn(k,k_t(t,1));
    c_t(t,1)=policyc1(index_k);
    l_t(t,1)=policyl1(index_k);
    k_t(t+1,1)=policyk1(index_k);
end

% compute implementability constraint
bla_t=zeros(T,1);
for t=1:T
    if sigma==1
        bla_t(t,1)=beta.^(t-1).*(theta-(1-theta).*l_t(t,1)./(1-l_t(t,1)));
    else
        bla_t(t,1)=beta.^(t-1).*(theta-(1-theta).*l_t(t,1)./(1-l_t(t,1))).*(c_t(t,1).^theta.*(1-l_t(t,1)).^(1-theta)).^(1-sigma);
    end
end

if sigma==1
    imc=sum(bla_t)-theta./c_t(1,1).*k(1,1);
else
    imc=sum(bla_t)-theta.*(c_t(1,1).^theta.*(1-l_t(1,1)).^(1-theta)).^(1-sigma)./c_t(1,1).*k(1,1);
end

imc

mu=.8*mu+.2*(mu-imc/2)
crit_mu=abs(imc)
crit_iter=crit_iter+1
end
    
% lambda
if sigma==1
    lambda_3=(exp((1-beta)*v0(index_k0))./(1-l0).^(1-theta)).^(1/theta)./c0-1;
else
    lambda_3=(((1-beta)*(1-sigma)*v0(index_k0)).^(1/(1-sigma))./(1-l0).^(1-theta)).^(1/theta)./c0-1;
end

lambda_3

figure(3); subplot(3,1,1); plot(k',v0'); title('value function from ramsey problem as function of k(kappa)'); grid on
figure(3); subplot(3,1,2); plot([1:50],[c_t(1:50)'; l_t(1:50)'; k_t(1:50)']); title('evolution of optimal c, l, k '); grid on
figure(3); subplot(3,1,3); plot([1:50],[kappa_t(1:50)'; tau_t(1:50)']); title('evolution of optimal kappa and tau'); grid on

toc