% This is FindingProblems5 - it is the same as FindingProblems4 except that
% it calls OF2Decisions5 which is different from OF2Decisions4.
% This is FindingProblems4 - see the Word document "Finding Problems" in the
% documentation folder. Version 3 differs from previous versions in order to save
% memory. Version 4 ignores the decisions in the final period which must
% all be zero and hence are not interesting.

% This version calls RollingDecsions5, OF1Decisions5, OF2Decisions5 which
% ignore the final period, Except for OF2Decisions5 these are all the same
% as the 4 versions.
% RollingDecsions3, OF1Decisions3 and OF2Decisions3
% differ from earlier versions in that they calculate things for just 1
% S value and just 1 r value respectively.
% It will produce five matrices, one for each of the risk aversion index r for 
% each of the four specifications and the rolling horizon S . The matrix entries will
% record the number of continuation decisions.
% These matrices will be denoted by ctxyz where x takes values r or S, y takes the values 1 or 2
% and z takes the values R (CRRA) or A (CARA). 
% We first need to read in the problems from problems.txt in this folder.
% CARE - in the problems.txt file in this directory all the p's are
% multiplied by 10 and all the rr's are multiplied by 100.
% I am today - because of lack on memory - returning just the number of
% continuation decisions rather than all of them. 
problemfile = fopen('problems.txt','r');
% Ignore first four lines
fgetl(problemfile);
fgetl(problemfile);
% read data
problems=single(cell2mat(textscan(problemfile,'%f %f %f %f %f %f')));
pall=problems(:,1);
hall=problems(:,2);
Tall=problems(:,3);
Lall=problems(:,4);
rrall=problems(:,5);
icfall=problems(:,6);
fclose(problemfile);
np=size(pall,1);
maxT=max(Tall);
Maxtotk=maxT*(maxT+1)/2;

nr=10;
rsr=1.0;  % CRRA
rdr=-0.1; % CRRA
rsa=0.0;  % CARA
rda=0.01; % CARA
fp=1;     % first problem
lp=np;     % last problem

% We have four specifications and two utility functions for the
% risk-aversion case. Let us do them one by one.

% First the Rolling Horizon
ctS=zeros(maxT-1,1);
pcctctS=zeros(np,1);
sdctS=zeros(np,1);
for pn=fp:1:lp
    p=pall(pn);
    h=hall(pn);
    T=Tall(pn);
    L=Lall(pn);
    rr=rrall(pn);
    icf=icfall(pn);
    for S=1:1:T-1
        % ctS is now the thing returned - the sum of the continuation decisions
        ctS(S) = RollingDecisions5(p,h,T,L,rr,icf,S);
    end
    pcctctS(pn)=100*mean(ctS)/(T*(T+1)/2);
    sdctS(pn)=std(ctS);
end

% Now the risk aversion cases. We need both OF1 and OF2.

% First OF1 and CRRA
ctr11=zeros(nr,1);
pcctctr11=zeros(np,1);
sdctr11=zeros(np,1);
for pn=fp:1:lp
    p=pall(pn);
    h=hall(pn);
    T=Tall(pn);
    L=Lall(pn);
    rr=rrall(pn);
    icf=icfall(pn);
    for j=1:1:nr
        r=rsr+(j-1)*rdr;
        % ctr11 is now the thing returned - the sum of the continuation
        % decisions
        ctr11(j) = OF1Decisions5(p,h,T,L,rr,icf,1,r);
    end
    pcctctr11(pn)=100*mean(ctr11)/(T*(T+1)/2);
    sdctr11(pn)=std(ctr11);
end

% Now OF1 and CARA
ctr12=zeros(nr,1);
pcctctr12=zeros(np,1);
sdctr12=zeros(np,1);
for pn=fp:1:lp
    p=pall(pn);
    h=hall(pn);
    T=Tall(pn);
    L=Lall(pn);
    rr=rrall(pn);
    icf=icfall(pn);
    for j=1:1:nr
         r=rsa+(j-1)*rda;
         % ctr12 is now the thing returned - the sum of the continuation
         % decisions
         ctr12(j) = OF1Decisions5(p,h,T,L,rr,icf,2,r);
    end
    pcctctr12(pn)=100*mean(ctr12)/(T*(T+1)/2);
    sdctr12(pn)=std(ctr12);
end

% Now OF2 and CRRA
ctr21=zeros(nr,1);
pcctctr21=zeros(np,1);
sdctr21=zeros(np,1);
for pn=fp:1:lp
    p=pall(pn);
    h=hall(pn);
    T=Tall(pn);
    L=Lall(pn);
    rr=rrall(pn);
    icf=icfall(pn);
    for j=1:1:nr
        r=rsr+(j-1)*rdr;
        % ctr21 is now the thing returned - the sum of the continuation
        % decisions
        ctr21(j) = OF2Decisions5(p,h,T,L,rr,icf,1,r);
    end
    pcctctr21(pn)=100*mean(ctr21)/(2^(T-1));
    sdctr21(pn)=std(ctr21);
end

% Now OF2 and CARA
ctr22=zeros(nr,1);
pcctctr22=zeros(np,1);
sdctr22=zeros(np,1);
for pn=fp:1:lp
    p=pall(pn);
    h=hall(pn);
    T=Tall(pn);
    L=Lall(pn);
    rr=rrall(pn);
    icf=icfall(pn);
    for j=1:1:nr
        r=rsa+(j-1)*rda;
        % ctr22 is now the thing returned - the sum of the continuation
        % decisions
        ctr22(j) = OF2Decisions5(p,h,T,L,rr,icf,2,r);
    end
    pcctctr22(pn)=100*mean(ctr22)/(2^(T-1));
    sdctr22(pn)=std(ctr22);
end

pcctall=[pcctctS,pcctctr11,pcctctr12,pcctctr21,pcctctr22];
sdall=[sdctS,sdctr11,sdctr12,sdctr21,sdctr22];
inandout=[pall,hall,Tall,Lall,rrall,icfall,pcctall,sdall];