tic clear all global tau k y_low theta_low theta_high gamma % parameters tau = 0.3; k = 4; y_low = 2; % truncate distribution at the top Fhigh = 1-1e-3; theta_low = theta_high = y_high = ygrid = y_low:(y_high-y_low)/1000:y_high; % gamma = 0.50; rr = 1; iter = 1; while (abs(rr)> 1e-3) && (iter < 1000); % v0 = fzero(@transversality,[0.2 1000]); [tt X] = ode45(@motion,[theta_low theta_high],[0 v0]); vv = X(:,2); eta = X(:,1); clear X % f = zeros(length(tt),1); for i=1:length(tt) f(i) = density(tt(i)); end % yy = cc = % r = zeros(1,length(tt)-1); for i=1:length(tt)-1 r(i) = (cc(i)-yy(i))*f(i)*(tt(i+1)-tt(i)); end rr=sum(r); % gamma = gamma + 0.001*sign(rr); iter = iter+1; end % plots figure(1); plot(yy,yy-cc); title('tax function (y-c as function of y)'); grid on figure(2); subplot(2,1,1); plot(yy,1-yy./(tt.^2)); title('marginal tax (as function of y)'); grid on subplot(2,1,2); plot(yy,(yy-cc)./yy); title('average tax (as function of y)'); grid on figure(3); subplot(2,1,1); plot(tt,eta); title('eta (as function of theta)'); grid on subplot(2,1,2); plot(tt,vv); title('v (as function of theta)'); grid on figure(4); subplot(2,1,1); plot(tt,cc); title('c (as function of theta)'); grid on subplot(2,1,2); plot(tt,yy); title('y (as function of theta)'); grid on figure(5); plot(ygrid,k*ygrid.^(-k-1)*y_low^k); title('density of income'); grid on