
% This program generates our benchmark model's saving rate and plots 
% its time path with the actual Japanese saving rate.
% See Figure 2 of Chen, Imrohoroglu, and Imrohoroglu in AER.

% The structure of the shooting algorithm is nicely summarized
% in Judd (1998), pp. 355-357. 

load datafile.txt ; % 45 by 9 data matrix, see below, starting with line 33

% parameters related to shooting algorithm
N = 2000;        % length of period used for forward shooting algorithm
epsilon = 0.001; % convegence criteria for shooting algorithm 
transp = 100;    % number of periods until convergence to the steady-state


% model parameters at the steady state

deltass = 0.10;         % depreciation rate for physical capital at the steady state;
beta    = 0.96327;      % utility discount factor     
theta   = 0.363;        % share of capital in production function
tauss   = 0.35;         % capital income tax rate at the steady state
govstar = 0.14;         % share of government spending in GNP at the steady state;
labor   = 31.82;        % exogenous labor input average over this time period
gammas  = 1.0315;       % steady state TFP factor growth rate (corresponds to 2% TFP growth rate)
popstar = 1.011976;     % steady state value of growth rate of time endowment
alpha   = 1.45;         % parameter for disutility of labor
T1      = 105;

b=1;                % starting year of simulation (b=1 refer to 1956)

% Loading exogenous time series from the .dat file
year    = datafile(b:45,1); 
gamma   = datafile(b:45,2);    % growth rate of TFP factor
pop     = datafile(b:45,3);    % population growth rate
gov     = datafile(b:45,4);    % government spending share
hh      = datafile(b:45,5);    % hours input
ee      = datafile(b:45,6);    % employment rate
delta   = datafile(b:45,7);    % depreciation rate
tau     = datafile(b:45,8);    % capital income tax
jsrate  = datafile(b:45,9);    % actual Japanese net national saving rate


% When feeding in steady state numbers use these:
%he      = labor*ones(45,1);
%pop     = popstar*ones(45,1);
%delta   = deltass*ones(45,1);
%tau     = tauss*ones(45,1);
%gov     = govstar*ones(45,1);
%gamma   = gammas*ones(45,1)

datlong  = length(year);


% Generating exogenous values for year beyond 2001
for j=datlong+1:N
    
    year(j)      = year(j-1)+1;
    gamma(j)     = gammas;
    tau(j)       = tauss;   
    delta(j)     = deltass; 
    pop(j)       = popstar; 
    gov(j)       = govstar; 
    end

    % Steady state values
    
    kh_star=((gammas-beta+deltass*beta*(1-tauss))/(beta*(1-tauss)*theta))^(1/(theta-1));

    A=(1-govstar)*kh_star^theta-kh_star*(gammas*popstar+deltass-1);

    hstar=T1/(1+alpha*A/((1-theta)*kh_star^theta));
    
    %pause

    cstar =A*hstar;


    kstar = kh_star*hstar;

    ystar = kstar^theta*hstar^(1-theta);

    save_ss = kstar*(gammas*popstar-1)/(ystar-deltass*kstar);  
    % steady state value of saving rate

    ky_star = kstar/ystar;    
    % steady state value of capital-output ratio
    
    % Calculation of transition path, using forward shooting algorithm
    c_L = cstar*0.0;            % initial guess for lower bound for consumption
    c_H = cstar*1.2;            % initial guess for upper bound for consumption
    % Follow Judd (p.357) and begin a shooting algorithm 

    % Setting initial capital stock
    ky_ini  =1.37;
    % k/y ratio at 1956           
    k0      = ky_ini^(1/(1-theta))*45.5;  % inital capital stock so that capital output ratio equal to data in 1956
    
   
    y       = 0;                         % dummy variable for convergence
    lc      = 0;                         % dummy for time

    while (y~=3 & lc < 100);
    
        y       =0;
    
        % initialize sequences to be calculated to zero
        kk      = zeros(N,1); 
        cc      = zeros(N,1);
        kkgg    = zeros(N,1);
        ccgg    = zeros(N,1);
        he      = zeros(N,1);
    
        kk(1,1) = k0;
    
        t       = 1;
        lc      = lc+1;
    
        % update guess for consumption at initial period
        cc(1,1) = (c_L+c_H)/2;

        c0      = cc(1,1);
    
    
       % disp([c_L c0 c_H]);
    
        % initial guess for h0
        X1=50.7;
        OPTIONS=optimset('display','off');
    
        % solve for h0, given k0 and c0
        [h0, fval, exitflag]=fsolve('fvalue_ini2', X1, OPTIONS, c0, alpha, theta, k0, T1);
    
    
        he(1,1) = h0; 
    

        % start shooting
        while t<N;
        
            t=t+1;
            kk(t,1) = ((1-delta(t-1,1))*kk(t-1,1)+(1-gov(t-1,1))*(kk(t-1,1)^theta*he(t-1,1)^(1-theta))...
            -cc(t-1,1))/(gamma(t-1,1)*pop(t-1,1));
        
            if kk(t,1) < 0;
               % disp('c_H will be lowered because kk is negative'), disp(t)
                y=1;
                %pause
                break;
            elseif imag(kk(t,1))~=0;
            %    disp('c_H will be lowered because kk is imaginary'), disp(t)
                y=1;
                break;
            end;
            
        
            % initial guess for c2
            X0= [cc(t-1,1), he(t-1,1)];
        
            % solve for optimal c and h for tomorrow 
            OPTIONS=optimset('display','off');
        
            [X2, fval, exitflag]=fsolve('fvalue2', X0, OPTIONS, cc(t-1,1), gamma(t-1,1), beta, tau(t,1), theta, delta(t,1), alpha, kk(t,1),T1);    

            if exitflag <1
                display('error')
                pause
            end

            %hist(t-1)=exitflag;
        
            cc(t,1) = X2(1);

            he(t,1) = X2(2);
        
            cgrow = cc(t,1)/cc(t-1,1);
            kgrow = kk(t,1)/kk(t-1,1);
            ccgg(t-1,1) = cgrow;
            kkgg(t-1,1) = kgrow;
        
            %if t>transp
            %pause
            %end

            if cgrow <=1 & t>transp;
                y=2;
           %     disp('consumption growth is too small'), disp(y), display(t);
                %pause
                break;
            elseif kgrow <=1 & t>transp;
                y=1;
            %    disp('capital growth is too small'), disp(y), display(t) ;
                %pause 
                break;
            end;
        
            %pause
        
        end;    % while t<N

        %pause
        
        % check convergence and update c_L (or c_H)
        if abs(cc(t,1)-cstar) < epsilon & abs(kk(t,1)-kstar) < epsilon*10 ;
            y=3;
            break;
        elseif y==2;
            c_L = cc(1,1); 
      %      disp('c_L increased');
        elseif y==1; 
            c_H = cc(1,1);
       %     disp('c_H lowered');
        else
      %      disp('number of iterations exceeded without convergence N and t are'), disp([N, t]);
            break;
    
        end;
    
        CHL = [c_L; c0; c_H]; 
    
    end; %(while y~=3..)
    
% calculating values in the model
N=max(find(kk));
kk=kk(1:N);
he=he(1:N);
%kk=kk(find(kk));
%he=he(find(he));
delta=delta(1:N);
cc=cc(1:N);
gov=gov(1:N);
tau=tau(1:N);


yy      = kk.^theta.*he.^(1-theta);                                  % GNP
rr      = (yy*theta./kk)-delta;                                      % rate of return for capital
rr      = rr.*(ones(N,1)-tau);                                       % after tax rate of return
savrate = (yy.*(ones(N,1)-gov)-cc-kk.*delta)./(yy-kk.*delta);        % saving rate
ky      = kk./yy;                                                    % capital output ratio


% save results
year    = [1956:2020];
result1 = [year', kk(1:65,1), cc(1:65,1), yy(1:65,1),ky(1:65,1),rr(1:65,1), he(1:65,1), savrate(1:65,1)];

plot(year(1:45),[jsrate savrate(1:45)]);
title('Figure 2. Benchmark Economy');
ylabel('Saving Rate');


% If you want to write the numerical results into a file, 
% remove the % signs from below:

%fid=fopen('result_inf.out', 'w');
%fprintf(fid,'%4i %8.4f  %8.4f %8.4f %7.4f %7.4f %7.4f %7.4f\n',result1');
%fprintf(fid,'\n');
%fprintf(fid,'\n');
%fclose (fid);


