%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%% % 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