%By David Love, Yale, using original code by Hugo Benitez-Silva, Summer 2001.
%Modified by Hugo Benitez-Silva, September 2001-March 2009.

%This program has the setup and the main program in one set of code. setupandbicrra1.m

%Setup contains some basic definitions and parameter values for the lifecycle problem

tt=60;                      %Number of finite time periods
periods=tt+1;
simu=500;
global vf i wgrid mu_r sig_r ar_r qa qw lnmu lnsig lnmx u1 ddelta bq gama beta intmeth wmax lntden unquad rp periods itt sut;
beta=.98;                   %Discount factor
bq=.1;                      %Bequest factor
gama=1.5;                   %CRRA parameter
ar_r=0;                     %Autoregressive parameter for lognormal AR(1)
lnparametersh=@lnparameters;
[mu_r,sig_r]=feval(lnparametersh,(1/beta),.3);
mu_r=mu_r*(1-ar_r);
sig_r=sig_r*(1-ar_r^2);

wn=100;                     %Number of grid point on w axis (0 is the (log(rmax)-rm)/rs1st grid point).
nrint=50;                   %Number of sample points used in expected value calculations.
nqint=10;                   %Number of quadrature abscissa for expected value.
wsp=1.3;                    %Determines degree of nonuniformity in spacing of w axis.
                            %grid points: wsp=1 gives uniformly spaced points, wsp>1.
                            %gives closer spacing near 0 and wider spacing as q increases.
rmax=5;                     %Upper bound of truncated lognormal density for returns.
rmin=0;                     %Lower bound of truncated lognormal density for returns.
macheps=1e-16;              %Tolerance for equality comparisons.
wmax=200000;                 %Upper bound for wealth.
wmin=1e-5;                  %Lower bound for wealth
told=1e-16;                 %Tolerance for equality comparisons.
btol=1e-2;                  %Tolerance for equality comparisons.feval(dufunc1h,cnew)+beta*(1-sut(periods-itt))*feval(dev_qw_myh,w,cnew)-sut(periods-itt)*bq*feval(dufunc1h,a)
ddelta=1e-7;                %Difference qoutient for numerical derivatives using superaccurate
                            %central differences.  Optimal values taken from table on p. 186
                            %from Rice, 1983, Numerical Methods, Software and Analysis, McGraw Hill.
if (wmin<ddelta);
    wmin=ddelta;
end;
intmeth=3;                  %0= monte carlo/low discrepancy integration.
                            %1= probability weighted integration.
                            %2= true expected value function.
                            %3= Gaussian quadrature.
unquad=1;                   %1= Do quadrature on interval (0,1) using probability integral
                            % transform method for untruncated lognormal.
                            %0= Do quadrature on interval(0,rmax) with truncated lognormal
                            %density.
qgaush=@qgaus;
if intmeth==3;
    if unquad==0;
        [qa,qw]=feval(qgaush,nqint,0,rmax);
    else;
        [qa,qw]=feval(qgaush,nqint,0,1);
    end;
end;

%Numbers and vectors to be used in ev_qw and other functions;
lntmu_rh=@lntmu_r;
lntsig_rh=@lntsig_r;
lntden_rh=@lntden_r;
rm=feval(lntmu_rh,rmax);
rs=feval(lntsig_rh,rmax);
ru=norm1((log(rmax)-rm)/rs);
rp=exp(rm+rs*norm2(ru*qa));
lnmx=rmax;
lntmu_rh=@lntmu_r;
lntsig_rh=@lntsig_r;
lnmu=feval(lntmu_rh,rmax);
lnsig=feval(lntsig_rh,rmax);
lntden=feval(lntden_rh,qa);
ppi=0;                      %Switch for parametric policy iteration.  This should be zero unless
                            %you are running ppi.gpr, but that program will automatically switch
                            %ppi=1.
%home=path(path,'/home...') %Path to directory to hold temporary working arrays for "virtual memory"
                            %versions of the programs.
%vmdir='/space1/jrust';
%readrows=50;                %Number of rows read in at each pass in virtual memory version.
if intmeth<3;
    wgrid=zeros(wn,1);
    wgrid(1)=wmin;
    i=2;
    while i<=wn;
        wgrid(i)=wgrid(i-1)+(wmax-wgrid(i-1))/(wn-i+1)^wsp;
        i=i+1;
    end
    wlb=[0;(wgrid(2:wn)+wgrid(1:wn-1))/2];
else;
    wgrid=zeros(5,1);
    wsp=1.5;
    wgrid(1)=wmin;
    i=2;
    while i<=size(wgrid,1);
        wgrid(i)=wgrid(i-1)+(wmax-wgrid(i-1))/(size(wgrid,1)-i+1)^wsp;
        i=i+1;
    end;
    divpt=wgrid;
    dwn=[floor(wn/2);floor(wn/4);floor(wn/6)];
    dwn=[dwn;wn-sum(dwn)];
    [wgrid,wts]=feval(qgaush,dwn(1),wmin,divpt(2));
    [wgrid0,wts0]=feval(qgaush,dwn(2),divpt(2),divpt(3));
    [wgrid1,wts1]=feval(qgaush,dwn(3),divpt(3),divpt(4));
    [wgrid2,wts2]=feval(qgaush,dwn(4),divpt(4),divpt(5));
    wgrid=[wgrid;wgrid0;wgrid1;wgrid2];
    wts=[wts;wts0;wts1;wts2];
end;

u1=[0.000996
0.000991
0.001001
0.001033
0.001081
0.001136
0.001192
0.001256
0.001327
0.001404
0.001485
0.001572
0.001670
0.001784
0.001915
0.002060
0.002216
0.002387
0.002573
0.002776
0.003003
0.003254
0.003522
0.003805
0.004107
0.004444
0.004825
0.005244
0.005708
0.006234
0.006845
0.007548
0.008327
0.009161
0.010046
0.011011
0.012078
0.013214
0.014394
0.015604
0.016786
0.018024
0.019482
0.021274
0.023376
0.025652
0.027986
0.030432
0.032968
0.035629
0.038427
0.041468
0.044936
0.049036
0.053854
0.059383
0.065545
0.072407
0.079904
0.088120];

sut=zeros(periods-1,1);
i=1;
while i<=(periods-1);
    sut(i)=u1(i);
    i=i+1;
end;
wealth=10000;

%include -a,cfunc in user library.


%bisectnew1h=@bisectnew1;
%evct1h=@evct1;
%devct1h=@devct1;
dufunc1h=@dufunc1;
%evc1h=@evc1;
%devc1h=@devc1;
%dev_qw_myh=@dev_qw_my;
%ev_qw_myh=@ev_qw_my;

%bicrra1.m: Program to compute the solution of the life-cycle
%consumption/savings problem via backward induction. CRRA utility.
%By David Love, Yale, using original code by Hugo Benitez-Silva, Summer 2001.
%Modified by Hugo Benitez-Silva, September 2001-September 2005.

t0=cputime;
format short;
%load setup.m;

%if intmeth==0;
%    disp('Program uses Monte Carlo/Low discrepancy integration');
%elseif intmeth==1;
%    disp('Program uses quadrature integration');
%end;

vf=zeros(wn,1);
vf1=vf;
c1=vf;
itt=1;
tol=1e-2;

i=1;
%bisectnew1h=@bisectnew1;
%bellman1tch=@bellman1tc;
%bellman1tch=@bellman2tc;
%bellmallch=@bellmallc;
devct2h=@devct2;
evct2h=@evct2;

while i<=size(wgrid,1);
   disp(i);
   %[c1(i),vf1(i)]=feval(bellman1tch,wgrid(i));
   c1(i)=fzero(devct2h,[ddelta wgrid(i)-(ddelta)]);
   vf1(i)=feval(evct2h,c1(i));
    
    i=i+1;
end;
t1=cputime;
disp('Last period Done in the following minutes')
disp((t1-t0)/60); 

vf=vf1;
vfcrnov10c=vf1;
c=c1;
cscrnov10c=c1;


%evc3h=@evc3;
dev_qw_myh=@dev_qw_my1;
ev_qw_myh=@ev_qw_my1;
%i=1;     
%fff(1)=feval(ev_qw_myh,5+ddelta);
%fff(2)=feval(ev_qw_myh,5-ddelta);
%vf1(2)=feval(dev_qw_myh,5);
%vf1(i+1)=feval(devc3h,wgrid(1)-ddelta);

devc3h=@devc3;
evc3h=@evc3;
while itt<=tt;
    i=1;
    disp(itt);
    while i<=size(wgrid,1);
        %[c1(i),vf1(i)]=feval(bellmallch,wgrid(i));
        c1(i)=fzero(devc3h,[ddelta wgrid(i)-(ddelta)]);
        vf1(i)=feval(evc3h,c1(i));
        i=i+1;
    end;
    vfcrnov10c=[vfcrnov10c;vf1];
    cscrnov10c=[cscrnov10c;c1];
    vf=vf1;
    vf1=zeros(wn,1);
    itt=itt+1;
end;

vfcr1=vfcrnov10c;
cscr1=cscrnov10c;
save bicrra.mat vfcr1 cscr1;
t1=cputime;
disp('Program Ends');
disp((t1-t0)/60);
disp('minutes');