% Carroll, Overland, and Weil (1997, 2000)
% Matlab programs 2003-2007 © Patrick Toche

function cow1
clear all; delete thediary.out; diary thediary.out;
disp('Patrick Toche©');
disp('The Carroll, Overland, and Weil (2000) model');
disp('The habit formation mechanism is defined in terms of consumption');
global rho theta delta;
rho=0.2; theta=0.05; delta=0.05;
fprintf('rho = %7.3f, theta = %7.3f, delta = %7.3f\n',rho,theta,delta);

%%% calibrating the initial steady state
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
global sigma psi gam Delta;
%Delta=3; sigma=5; psi=sigma-Delta; gam=psi/(sigma-1);
Delta=3; psi=4; sigma=psi+Delta; gam=psi/(sigma-1);
%sigma=9; gam=0.75; psi=gam*(sigma-1); Delta=sigma-psi;
%sigma=3; psi=1.00; gam=psi/(sigma-1); Delta=sigma-psi;
fprintf('sigma = %7.3f, psi = %7.3f, gamma = %7.3f, Delta = %7.3f\n',sigma,psi,gam,Delta);
global g A;
g=0.02; A=delta+theta+Delta*g;
global xss yss y0 s dsdg;
s=(g+delta)/A;
dsdg=(delta+theta-delta*Delta)/(delta+theta+g*Delta)^2;
fprintf('g = %7.3f, A = %7.3f, s = %7.3f, dsdg = %7.3f\n',g,A,s,dsdg);
xss=1+(A-delta-theta)/(rho*Delta);
yss=xss/(A-delta-rho*(xss-1));
%y0=1.5*yss;
y0=0.5*yss;
fprintf('xss = %7.3f, yss = %7.3f, y0 = %7.3f\n',xss,yss,y0);

%%% checking that the computed steady state annulls the ODE equations
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%ess=zeros(2,1); ess(1)=xss; ess(2)=yss; odess=feval(@ODE,1,ess); fprintf('%5.6f\n',odess);
%jss=zeros(2,2); jss=feval(@J,1,ess); tr=trace(jss); dt=det(jss); fprintf('tr = %5.6f, det = %5.6f\n',tr,dt);

%%% computing the dynamics using the bvp4c routine
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
global T0 T N NMax RelTol AbsTol guess;
T0=10; T=input('Horizon Length:  '); Tmax=min(50,T);
N=10; NMax=500; RelTol=1e-4; AbsTol=1e-2;
guess=[xss y0];
%options = bvpset('Vectorized','on','Stats','on','RelTol',RelTol,'AbsTol',AbsTol','NMax',NMax);
options = bvpset('FJacobian',@J,'Vectorized','on','Stats','on','RelTol',RelTol,'AbsTol',AbsTol','NMax',NMax);
solinit = bvpinit(linspace(0,T0,N),guess);
sol = bvp4c(@ODE,@BC,solinit,options);
t = linspace(0,T0);
e = deval(sol,t);
e(3,:) = (A-delta-theta+psi*rho.*(e(1,:)-1))/sigma;  %growth rate of consumption
e(4,:) = A-delta-e(1,:)./e(2,:); %growth rate of capital
e(5,:) = 1-e(1,:)./(A.*e(2,:)); %saving rate
fprintf('With T = %g, x(0) = %7.7f, y(0) = %7.7f\n',T0,e(1,1),e(2,1));
fprintf('With T = %g, x(T) = %7.7f, y(T) = %7.7f\n',T0,e(1,end),e(2,end));

%%% drawing some figures
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

figure(101);
clf;
plot(t,e(1,:),'-',t(end),e(1,end),'o');
title('consumption/habit ratio');
xlabel('t');
ylabel('x(t)');
hold on;
drawnow;

figure(102);
clf;
plot(t,e(2,:),'-',t(end),e(2,end),'o');
title('capital/habit ratio');
xlabel('t');
ylabel('y(t)');
hold on;
drawnow;

%%% use continuation to adjust end-point
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
for TT = T0+1:T
solinit = bvpinit(sol,[0 TT]);
sol = bvp4c(@ODE,@BC,solinit,options);
t = linspace(0,TT);
e = deval(sol,t);
e(3,:) = (A-delta-theta+psi*rho.*(e(1,:)-1))/sigma;  %growth rate of consumption
e(4,:) = A-delta-e(1,:)./e(2,:); %growth rate of capital
e(5,:) = 1-e(1,:)./(A.*e(2,:)); %saving rate
fprintf('With T = %g, x(0) = %7.7f, y(0) = %7.7f\n',TT,e(1,1),e(2,1));
fprintf('With T = %g, x(T) = %7.7f, y(T) = %7.7f\n',TT,e(1,end),e(2,end));
fprintf('With T = %g, s(0) = %7.7f, s(T) = %7.7f\n',TT,e(5,1),e(5,end));

%%% completing the figures
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
figure(101);
plot(t,e(1,:),'-',t(end),e(1,end),'o');
figure(102);
plot(t,e(2,:),'-',t(end),e(2,end),'o');
drawnow;
end
hold off;

%%% drawing more figures
%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%figure(1)

figure(1);
clf;
plot(t,e(1,:),'-',t(end),e(1,end),'o',t(1),e(1,1),'o');
grid on;
title('consumption/habit ratio');
xlabel('t');
ylabel('x(t)');
hold on;
drawnow;

figure(2);
clf;
plot(t,e(2,:),'-',t(end),e(2,end),'o',t(1),y0,'o',t(1),e(2,1),'o');
grid on;
title('capital/habit ratio');
xlabel('t');
ylabel('y(t)');
hold on;
drawnow;

figure(3);
clf;
plot(t,e(3,:),'-',t(1),e(3,1),'o',t(end),e(3,end),'o');
grid on;
title('consumption growth');
xlabel('t');
ylabel('gc');
hold on;
drawnow;

figure(4);
clf;
plot(t,e(4,:),'-',t(1),e(4,1),'o',t(end),e(4,end),'o');
grid on;
title('capital growth');
xlabel('t');
ylabel('gk');
hold on;
drawnow;

figure(5);
clf;
plot(t,e(5,:),'-',t(1),e(5,1),'o',t(end),e(5,end),'o');
grid on;
title('saving rate');
xlabel('t');
ylabel('s/y');
hold on;
drawnow;

figure(6);
clf;
plot(e(2,:),e(1,:),'-',e(2,1),e(1,1),'o',e(2,end),e(1,end),'o');
grid on;
title('consumption/habit against capital/habit ratios');
xlabel('k/h');
ylabel('c/h');
hold on;
drawnow;

figure(7);
clf;
plot(e(2,:),e(3,:),'-',e(2,1),e(3,1),'o',e(2,end),e(3,end),'o');
grid on;
title('consumption growth against capital/habit ratios');
xlabel('k/h');
ylabel('gc');
hold on;
drawnow;

figure(8);
clf;
plot(e(2,:),e(4,:),'-',e(2,1),e(4,1),'o',e(2,end),e(4,end),'o');
grid on;
title('capital growth against capital/habit ratios');
xlabel('k/h');
ylabel('gk');
hold on;
drawnow;

figure(9);
clf;
plot(e(2,:),e(5,:),'-',e(2,1),e(5,1),'o',e(2,end),e(5,end),'o');
grid on;
title('saving rate against capital/habit ratios');
xlabel('k/h');
ylabel('s/y');
hold on;
drawnow;

figure(10);
clf;
plot(e(3,:),e(5,:),'-',e(3,1),e(5,1),'o',e(3,end),e(5,end),'o');
grid on;
title('growth rate of consumption against saving rate');
xlabel('gc');
ylabel('s/y');
hold on;
drawnow;

hold off; diary off; save('cow.mat');

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
function rr = ODE(t,e)
global sigma psi rho theta A delta;
rr = [((A-delta-theta)-(sigma-psi)*rho.*(e(1,:)-1)).*e(1,:)/sigma;
      (A-delta-rho*(e(1,:)-1)).*e(2,:)-e(1,:)];

function rr = BC(e0,eT)
global  y0 yss;
rr = zeros(2,1);
rr(1) = eT(2)-yss;
rr(2) = e0(2)-y0;

function j = J(t,e)
global sigma psi rho theta A delta xss yss;
j = zeros(2,2);
j(1,1) = -(sigma-psi)*rho*xss/sigma;
j(1,2) = 0;
j(2,1) = -1;
j(2,2) = xss/yss;

function rr = U(c,z)  %utility function
global sigma psi;
rr = c^(1-sigma)*z^(psi);