% This is estimation 11. I am going to 'subscript' the variable od.
% od1r for OF=1 and rt=1
% od1a for OF=1 and rt=2
% od2r for OF=2 and rt=1
% od2a for OF=2 and rt=2
% ods  for OF=5

% This is estimation10. It re-does the graphs.
% Estimation 9 got rid of the stupid bug in all the previous versions - 
% by inserting the following lines (now lines 61 to 63):
% if lqv(k)>ctv(k)
%           EV(k)=lqv(k);
%        end.
% It now produces the same answers for a risk-neutral subject on both OF1
% and OF2.

% estimation5 differs from estimation4 in that it does OF=2.

% estimation4 differs from estimation3 in that it does both CRRA and CARA.

tic
% There are two different objective functions specified by the variable OFI:
% OFI=1 for sum of utilities.
% OFI=2 for utility of sum. 
% OFI=3 for both the above
% OFI=4 for rolling horizon.
% OFI=5 for all three.
% CARE OFI is telling it what to do. OF is specifying the Objective
% Function.
OFI=5;

% There are two different ways to find the solution:
% If EM=1 graphs the percentage of correct decisions,
% If EM=2 it finds the maximum likelihood estimate,
% If we want it to do both put EM=3.
EM=1;

% There are two different possible functional forms:
% rt=1 is CRRA
% rt=2 is CARA
% This version does both.
% I am using it for the main analysis 
% I am considering only problem set 1 - in Problems1.txt
% In each of which T=15 and there are 15 problems.
% This is my first bash at trying to find the best risk-attitude.
% Nodes are numbered sequentially from the start
% In the problems1.txt file the variables are as follows:
% 1: p (the probability of jumping up)
% 2: h (the jump size)
% 3: T (the total number of periods before compulsory liquidation)
% 4: minimum time
% 5: maximum time
% 6: L (the liquidation value)
% 7: rr  (the rate of return)
% 8: icf (the initial cash flow)
% What I think we need to produce, for each problem in the relevant input
% file, is the optimal decisions at each node in the tree for any given
% risk attitude. This will be a set of 1's and 0's for all (T+1)*T/2 nodes.
% I should write a function for this.

fsn=66;               % first subject to be analysed
lsn=66;               % last subject to be analysed
ns=74;               % number of subjects in main


restimate=zeros(ns,2);
maxcorrect=zeros(ns,2);


% I am just going to calculate 10 different risk indices
rsr=1.0;                   % this is the first CRRA risk aversion parameter plotted
rdr=-0.01;                 % this is the increment in the CRRA risk aversion parameter
rsa=0.0;                   % this is the first CARA risk aversion parameter plotted
rda=0.0005;                 % this is the increment in the CARA risk aversion parameter
nr=90;                     % this is the number of risk aversion parameters considered

% I first read in the problems.
% CARE - in the problems1.txt file in this directory all the p's are
% multiplied by 10 and all the rr's are multiplied by 100.
problemfile = fopen('problems1.txt','r');
% Ignore first three lines
fgetl(problemfile);
fgetl(problemfile);
fgetl(problemfile);
% read data
y=textscan(problemfile,'%d %d %d %d %d %d %d %d');
problems=single(cell2mat(y));
p=0.1*problems(:,1);
h=problems(:,2);
T=problems(:,3);
TT=T(1); % care we may have to change this if we are going to run sessions with different Ts.
L=problems(:,6);
rr=0.01*problems(:,7);
icf=problems(:,8);
fclose(problemfile);
np=size(p,1);
tott=TT;                   % these are the total number of t nodes
totk=TT*(TT+1)/2;          % these are the total number of k nodes
totl=2^TT-1;               % these are the total number of l nodes
totm=2^(TT-1)-1;           % these are the total number of l nodes in the penultimate period
% OK - we have read in the problems

% now read in the data.
nd=zeros(ns,1);    % this is the number of decisions by subjects.
%reading in the data of all the subjects.
load('DecsData');
sn=DecsData(:,1); 
pn=DecsData(:,2);
period=DecsData(:,3); % the period of the specific problem
kn=DecsData(:,4);  % this notation refers to the k node reached by the subjects (for OF1)
               % I will use ln for the l node reached by the subject.
ln=DecsData(:,5); % this notation refers to the l node reached by the subjects (for OF2)
ad=DecsData(:,6);
nl=size(sn,1); % this refers to the total number of nodes reached by all the subjects.

if ((OFI==1) || (OFI==2) || (OFI==3) || (OFI==5));
    if OFI==1
        OFF=1;
        OFL=1;
    end
    if OFI==2
        OFF=2;
        OFL=2;
    end
     if OFI==3
        OFF=1;
        OFL=2;
     end
    if OFI==5
        OFF=1;
        OFL=2;
    end
    for OF=OFF:1:OFL
        for rt=1:1:2
            if rt==1
                rs=rsr;
                rd=rdr;
            end
            if rt==2
                rs=rsa;
                rd=rda;
            end
      
%now we need to find the t value corresponding to any k value
kk=0;
for  tt=1:1:tott
    for k=1:1:tt
        kk=kk+1;
        ckt(kk)=tt;
    end
end

% Now generate od for OF=1.
if ((OFI==1) || (OFI==3) || (OFI==5))
    if rt==1
        od1r=-9*ones(totk,np*nr);
    end
    if rt==2
        od1a=-9*ones(totk,np*nr);
    end
    for rn=1:1:nr
        r=rs+(rn-1)*rd;
        for i=1:1:np
            pi=p(i);
            hi=h(i);
            Ti=T(i);
            Li=L(i);
            rri=rr(i);
            icfi=icf(i);
            cn=(rn-1)*np+i;
            if rt==1
                [od1r(:,cn)] = OptimalDecisions(rt,r,pi,hi,Ti,Li,rri,icfi,1);
            end
            if rt==2
                [od1a(:,cn)] = OptimalDecisions(rt,r,pi,hi,Ti,Li,rri,icfi,1);
            end
        end
    end
end

% now generate od for OF=2 This will have 15 times the number of risk indices columns 
%(one for each problem for each risk index)
% and 2^(T-1)-1 rows (one for each of the l 'nodes').
if ((OFI==2) || (OFI==3) || (OFI==5))
    if rt==1
        od2r=-9*ones(totl,np*nr);
    end
    if rt==2
        od2a=-9*ones(totl,np*nr);
    end
    % we need first to find the upl and dnl values.
    upl=zeros(totm,1);
    dnl=zeros(totm,1);
    %[od,upl,dnl,EV,lqv,ctv]=OptimalDecisions(2,0.000,0.5,1,6,15,1.1,0,2);
    %[od,upl,dnl,EV,lqv,ctv]=OptimalDecisions(2,0.1,0.5,1,6,15,1.1,0,2);
    [od,upl,dnl,acf,EV,lqv,ctv]=OptimalDecisions(1,1.0,0.5,1,15,10,1.1,0,2);
    
    % now I need to find the l node reached by the subject
    ln=zeros(nl,1);
    for l=1:1:nl
        if kn(l)==1
            ln(l)=1;
        end
        if kn(l)>1
            if kn(l)-kn(l-1)==ckt(kn(l-1))
                ln(l)=upl(ln(l-1));
            end
             if kn(l)-kn(l-1)==ckt(kn(l-1))+1
                ln(l)=dnl(ln(l-1));
            end
        end
    end
    ss=zeros(totl,nr);
    for rn=1:1:nr
        r=rs+(rn-1)*rd;
        for i=1:1:np
            pi=p(i);
            hi=h(i);
            Ti=T(i);
            Li=L(i);
            rri=rr(i);
            icfi=icf(i);
            cn=(rn-1)*np+i;
            % the following line calculates the optimal decisions for
            % particular problem and a particular risk index.
            if rt==1
                [od2r(:,cn)] = OptimalDecisions(rt,r,pi,hi,Ti,Li,rri,icfi,2);
            end
            if rt==2
                [od2a(:,cn)] = OptimalDecisions(rt,r,pi,hi,Ti,Li,rri,icfi,2);
            end
            %ss(:,rn)=ss(:,rn)+abs(od(:,cn)-od(:,i));
        end
        
    end
end



% Now the next bit only applies if the Estimation Method required is 1.
if ((EM==1) || (EM==3))
% now we start seeing how close the actual decisions are to the optimal for the
% various risk parameters
% I shall use the variable cd (correct decisions) to indicate this: for each
% subject and each risk parameter. So I will produce a matrix with ns (the number of subjects) 
% columns and nr rows.
% I need new notation. Not just one blanket cd and pcd for both OFs and for
% both rts. But let us not keep both cd and pcd. Let cd initially be the
% number correct and then make it the percentage correct. So cdxy denotes
% the number/percentage correct for OF x (1 or 2) and for rt y (r or a)
if ((OF==1) && (rt==1)) 
    cd1r=zeros(nr,ns);
end
if ((OF==1) && (rt==2))
    cd1a=zeros(nr,ns);
end
if ((OF==2) && (rt==1))
    cd2r=zeros(nr,ns);
end
if ((OF==2) && (rt==2))
    cd2a=zeros(nr,ns);
end

% We need to count the number of decisions of each subject.
% Let  us call this vector nd - with ns rows and 1 column.
nd=zeros(ns,1);
for l=1:1:nl
    nd(sn(l))=nd(sn(l))+1;
end
if OF==1
    for l=1:1:nl
        for rn=1:1:nr
            if rt==1
                if ad(l)==od1r(kn(l),(rn-1)*np+pn(l))
                    cd1r(rn,sn(l))=cd1r(rn,sn(l))+1;
                end
            end
             if rt==2
                if ad(l)==od1a(kn(l),(rn-1)*np+pn(l))
                    cd1a(rn,sn(l))=cd1a(rn,sn(l))+1;
                end
             end
        end
    end
end

if OF==2
    for l=1:1:nl
        for rn=1:1:nr
            if rt==1
                if ad(l)==od2r(ln(l),(rn-1)*np+pn(l)) 
                    cd2r(rn,sn(l))=cd2r(rn,sn(l))+1;
                end
            end
            if rt==2
                if ad(l)==od2a(ln(l),(rn-1)*np+pn(l)) 
                    cd2a(rn,sn(l))=cd2a(rn,sn(l))+1;
                end
            end
        end
    end
end

for snn=fsn:1:lsn
    for rn=1:1:nr
        if OF==1 && rt==1
            cd1r(rn,snn)=int8(100*cd1r(rn,snn)/nd(snn));
        end
        if OF==1 && rt==2
            cd1a(rn,snn)=int8(100*cd1a(rn,snn)/nd(snn));
        end
        if OF==2 && rt==1
            cd2r(rn,snn)=int8(100*cd2r(rn,snn)/nd(snn));
        end
        if OF==2 && rt==2
            cd2a(rn,snn)=int8(100*cd2a(rn,snn)/nd(snn));
        end
    end
end
end


% Now the next bit is if the Estimation Method required is 2. 
if ((EM==2) || (EM==3))
    if rt==1
        lb=rsr+(nr-1)*rdr;
        ub=rsr;
        r0=(lb+ub)/2;
    end
    if rt==2
        lb=rsa;
        ub=rsa+(nr-1)*rda;
        r0=(lb+ub)/2;
    end
options = psoptimset('MaxIter',1000);
    for snn=fsn:1:lsn   
        if OF==1
            [xval,fval,exitflag]=patternsearch(@(r) NegativeOfPercentageCorrect(r,rt,h,icf,L,p,rr,T,ad,nl,kn,pn,sn,snn,1),[r0],[],[],[],[],lb,ub,[],options);
        end
        if OF==2
            [xval,fval,exitflag]=patternsearch(@(r) NegativeOfPercentageCorrect(r,rt,h,icf,L,p,rr,T,ad,nl,ln,pn,sn,snn,2),[r0],[],[],[],[],lb,ub,[],options);
        end
        
        restimate(snn,rt+2*(OF-1))=xval;
        maxcorrect(snn,rt+2*(OF-1))=-fval;
    end
end
    end
    end
end

if ((OFI==4) || (OFI==5)) 
%now generate od for Rolling Horizon case
ods=-9*ones(totk,np*(TT-1));
for S=1:1:TT-1
    for i=1:1:np
    pi=p(i);
    hi=h(i);
    Ti=T(i);
    Li=L(i);
    rri=rr(i);
    icfi=icf(i);
    cn=(S-1)*np+i;
    [ods(:,cn),D]=RollingDecisions(S,pi,hi,Ti,Li,rri,icfi);
    end
end
% rod14 is the decisions for someone with the correct horizon
%rod14=od(:,196:210);
% now we count the number of decisions and the number corrrect
cds=zeros(TT-1,ns);
% It might also again be useful to count the number of decisions of each subject.
% Let  us call this vector nd - with ns rows and 1 column. As I am calling
% it again I should reset it to zero.
nd=zeros(ns,1);
for l=1:1:nl
    nd(sn(l))=nd(sn(l))+1;
end
for l=1:1:nl
    for S=1:1:TT-1
        if ad(l)==ods(kn(l),(S-1)*np+pn(l)) 
            cds(S,sn(l))=cds(S,sn(l))+1;
        end
    end
end

for snn=fsn:1:lsn
    for S=1:1:TT-1
        cds(S,snn)=int8(100*cds(S,snn)/nd(snn));
    end
end

if ((OFI==1) || (OFI==5))
    b=-13/nr;
    a=14-b;
    for rn=1:1:nr
        xval(rn)=a+b*rn;
    end
    
    for S=1:1:TT-1
        xvals(S)=S;
    end
    for snn=fsn:1:lsn
         figure('Name',['Percentage Correct Subject ',num2str(snn)],'NumberTitle','off')
         xlabel('horizon and risk attitude');
         ylabel('percentage correct');
         plot(xval,cd1r(1:nr,snn)','LineStyle','-','LineWidth',1);
         hold on
         plot(xval,cd1a(1:nr,snn)','LineStyle','--','LineWidth',3);
         hold on
         plot(xvals,cds(1:TT-1,snn)','LineStyle','-','LineWidth',2);
         axis([1 15 0 100])
         plot(xval,cd2r(1:nr,snn)','Marker','x','MarkerSize',2);
         hold on
         plot(xval,cd2a(1:nr,snn)','LineStyle','-.','LineWidth',1);
         hold on
         hold on
         saveas(gcf,['Percentage Correct Subject ',num2str(snn),'.jpg']);
         hold off
    end
         hold off
    end
end

save cd1r;
save cd1a;
save cd2r;
save cd2a;
save cds;
save nd;


% I do not know what the following does
% if ((EM==2) || (EM==3)) 
%     rncorrect=zeros(ns,4);
%     for snn=1:1:ns
%         rncorrect(snn,1) = -NegativeOfPercentageCorrect(1.0,1,h,icf,L,p,rr,T,ad,nl,kn,pn,sn,snn,1);
%         rncorrect(snn,2) = -NegativeOfPercentageCorrect(0.0,2,h,icf,L,p,rr,T,ad,nl,kn,pn,sn,snn,1);
%         rncorrect(snn,3) = -NegativeOfPercentageCorrect(1.0,1,h,icf,L,p,rr,T,ad,nl,ln,pn,sn,snn,2);
%         rncorrect(snn,4) = -NegativeOfPercentageCorrect(0.0,2,h,icf,L,p,rr,T,ad,nl,ln,pn,sn,snn,2);
%     end
% maxcorrect=[maxcorrect,rncorrect];
% end
toc


