close all; clear all; clc; tic;
% Instruction:
% Set SWT = 1; in line 35 or around it.
% Choose Computational and Deep parameters in lines 10-22.
% Some computation is unstable, for which error messages are surpressed.
% This is especially relevant to the case with eta < 1.0.


%% A: Parameters;
% Computational;
mxW     = 2.0;% 8.0;%       % maximum W;
stW     = 30;              % step size for W; stW & stR define resolution.
stR     = stW;              % step size for R;

% Deep;
rho     = 0.04;             % discount rate;
al2     = 0.0211;% 0.0117;% 0.02;% 0.050;% 0.061;% 0.019;%   % alpha~;
si2     = 0.0604;% 0.0112;% 0.20;% 0.040;% 0.042;% 0.060;%   % sigma~;
eta     = 5.0;% 0.7;%       % elasticity of substitution btw C and R;
phi     = 1.0;              % relative importance of R in pi;
kkk     = 0.0;              % marginal cost of conversion; 
                            % kkk>0 is not supported by this code.

% Non-Deep;
alp     = (1-1/eta)*(al2-si2*si2/2/eta);    % alpha;
sig     = (1-1/eta)*si2;                    % sigma;
sig2    = sig*sig;                          % sigma^2;
aa1     = 1/2 - alp/sig2;                   % a1 in beta;
aa2     = sqrt(aa1^2 + 2*rho/sig2);         % a2 in beta;
bet     = aa1 + aa2;                        % beta>1 +ve root;
ph2     = bet/(bet-1)*phi/rho/(1/(rho-alp)+kkk);                % phi2;
ph3     = ( phi/(rho*(bet-1)) )*ph2^(-bet);                     % phi3;
Bprim   = @(xxx)( xxx.^((bet-1)/eta) )./( (1-xxx).^(bet/eta) ); % B'(R)/ph3

SWT     = 1;                                % =1 in order not to save.
tmpN    =[al2*1000, si2*1000, eta*10];      % for saved file name.
KeyParam=[al2, si2, eta];                   %


%% B: Barrier;
www     = 0.0:(mxW/stW):mxW;                % vector of W;
zzz     = www.^(1-1/eta);                   % vector of Z;
rr1     = 1./( (zzz/ph2).^(  eta) +1 );     % optimal R
rr2     = 1./( (zzz/phi).^(  eta) +1 );     % myopic R
barr    = spline(www, rr1);                 % store barrier as a spline;

if SWT == 1;
figure('Name','Barrier on (R,W) plane');
temp = plot(rr1, www, rr2, www); axis([0,1,0,mxW]); grid on;
set(temp, 'LineWidth', 2); set(temp(2), 'LineStyle', '-.');
legend('stochastic case (eqn(  ))',  'myopic case     (eqn(  ))');
ylabel('Technology Level, W');  xlabel('Reserved Environment, R');
title(['Optimal vs. Myopic R: \eta=', num2str(eta)], 'FontSize', 14);
text(0.75, mxW*0.75,'Conversion Region');
text(0.02, mxW*0.25,'Non-Conversion Region');
else % if SWT == 0;
    tmpS=sprintf('Fig1Bar_A%02.0fS%02.0fE%02.0f', tmpN);
    save(tmpS, 'rr1',      'www',      'rr2',      'KeyParam');
end; % if SWT == 1;


%% C: Values;
%[C0] Create meshgrid, etc.
www         = 0.0:(mxW/stW):mxW;            % row vector;
rrr         = 0.0:( 1 /stR):1.0;            % col vector;
[RRR, WWW]  = meshgrid(rrr, www);           % matrices that copies rrr/www;
zzz         = www.^(1-1/eta);               % row vector;
ZZZ         = WWW.^(1-1/eta);               % matrix that copies rows; 

%[C1] B(R);
    bbr     = 0.0*ones(size(rrr));
for iiR = 1:1:numel(rrr)-1;                 % loop for integration
    warning('off', 'MATLAB:quad:MaxFcnCount'     );
    warning('off', 'MATLAB:quad:ImproperFcnValue'); % disp(num2str(iiR));
    bbr(iiR)= quad(Bprim, 0, rrr(iiR), 1e-12);      % adaptive simpson rule
    warning('on',  'MATLAB:quad:MaxFcnCount'     );
    warning('on',  'MATLAB:quad:ImproperFcnValue');
end; % for iiR = 2:1:numel(rrr)-1;
bbr         = ph3.*bbr;

%[C2] Option Value in NonConversion Region;
[BBR, WWW]  = meshgrid(bbr, www);           % B(R);
OpVal       = (  ZZZ.^bet  ).*(  BBR  );    % 3rd term: Option value;
tempOpVal   = OpVal;

for iiW = 1:1:numel(www);
    for iiR = 1:1:numel(rrr);
        if rrr(iiR) > ppval(www(iiW), barr);
            OpVal(iiW,iiR)  = NaN;          % for conversion region;
        end; 
    end; % for iiR = 1:1:numel(rrr);
end; % for iiW = 1:1:numel(www);

if SWT == 1;
figure('Name','Option Value in NonConversion Region ');
 meshc(RRR, WWW, OpVal); grid on; 
%contour(RRR, WWW, OpVal, 40); grid on; 
ylabel('Technology Level, W');  xlabel('Reserved Environment, R');
title(['Option Value in NonConversion Region: \eta=', num2str(eta)],...
       'FontSize', 14); 
else % if SWT == 1;
    tmpS=sprintf('Fig2OpV_A%02.0fS%02.0fE%02.0f', tmpN);
    save(tmpS, 'RRR',      'WWW',      'OpVal',      'KeyParam');
end; % if SWT == 1;
% toc; break;                                                       %%% <=
  
%[C3] NonOption Value in NonConversion Region; 
FFn     = ZZZ.*(1-RRR).^(1-1/eta)./(  (rho-alp)*(1-1/eta)  ) ...
        + phi.*(  RRR).^(1-1/eta)./(  (rho    )*(1-1/eta)  ); %
tempFFn = FFn;

tmpIdx  = NaN*ones(2,numel(www));
for iiW = 1:1:numel(www);
    tempSwtich  = 0;
    for iiR = 1:1:numel(rrr);
        if rrr(iiR) > ppval(www(iiW), barr);
            FFn(iiW,iiR)   = NaN;   % NaN for conversion region;
            if tempSwtich == 0;     % pick up the initial time only;
                tmpIdx(:,iiW)   = [iiW; iiR]; % recording cell position;
                tempSwtich      = 2.718281828;
            end; % if tempSwtich == 0;
        end; % if rrr(iiR) > ppval(www(iiW), Barr);
    end; % for iiR = 1:1:numel(rrr);
end; % for iiW = 1:1:numel(www);
if      eta < 1;        tmpIdx(:,1)     = [1; 0];
elseif  eta > 1;        tmpIdx(:,1)     = [1; 1];
elseif  eta ==1;        tmpIdx(:,1)     = [1; ph2/(1+phi2)];       
end; % this is to avoid inf or /0 for the first pair;

if SWT == 1;
figure('Name','NonOption Value in NonConversion Region');
meshc(RRR, WWW, FFn); grid on; 
ylabel('Technology Level, W');  xlabel('Reserved Environment, R');
title(['NonOption Value in NonConversion Region: \eta=', num2str(eta)],...
       'FontSize', 14); 
else % if SWT == 1;
    tmpS=sprintf('Fig3Nop_A%02.0fS%02.0fE%02.0f', tmpN);
    save(tmpS, 'RRR',      'WWW',      'FFn',      'KeyParam');
end; % if SWT == 1;

%[C4] Total Value in NonConversion Region;
FF0     = FFn + tempOpVal;
if SWT == 1;
figure('Name','Total Value in NonConversion Region');
meshc(RRR, WWW, FF0); grid on; 
ylabel('Technology Level, W');  xlabel('Reserved Environment, R');
title(['Total Value in NonConversion Region for \eta=', num2str(eta)], ...
       'FontSize', 14); 
else % if SWT == 1;
    tmpS=sprintf('Fig4NoC_A%02.0fS%02.0fE%02.0f', tmpN);
    save(tmpS, 'RRR',      'WWW',      'FF0',      'KeyParam');
end; % if SWT == 1;

%[C5] Total value in Conversion Region;
    br_     = 0.0*ones(size(rr1));
    warning('off', 'MATLAB:quad:MaxFcnCount'     );
    warning('off', 'MATLAB:quad:ImproperFcnValue'); % disp(num2str(iiR));
for iiW = 1:1:numel(rr1);               % loop for integration
    br_(iiW)= quad(Bprim, 0, rr1(iiW), 1e-12); % adaptive simpson's rule;
end; % for iiW = 2:1:numel(rr1);
    warning('on',  'MATLAB:quad:MaxFcnCount'     );
    warning('on',  'MATLAB:quad:ImproperFcnValue');
br_     = ph3*br_;                      % B(R) along the barrier;
opVa_   = (  zzz.^bet  ).*(  br_  );    % 3rd term along the barrier;

f0_     = zzz.*(1-rr1).^(1-1/eta)/(  (rho-alp)*(1-1/eta)  ) ...
        + rr1.^(1-1/eta)*phi/(   rho*(1-1/eta)  )           ...
        + opVa_;                        % total value along the barrier;
[~, F0_]= meshgrid(rrr, f0_);           % convert f0_ vec to mat(mesh);
[RRR, RR1]  = meshgrid(rrr, rr1);       % convert rr1 vec to mat(mesh);
cc      =  kkk/(1-1/eta)*ZZZ.*(  (1-RRR).^(1-1/eta) ... % 
                                -(1-RR1).^(1-1/eta)  ); % conversion cost;
FF1     = F0_ + cc;
tempFF1 = FF1;

for iiW = 1:1:numel(www);
    for iiR = 1:1:numel(rrr);
        if rrr(iiR) <= ppval(www(iiW), barr);
            FF1(iiW,iiR)    = NaN;      % NaN for conversion region;
            cc(iiW,iiR)     = NaN;
        end; % if rrr(iiR) <= ppval(www(iiW), barr);
    end; % for iiR = 1:1:numel(rrr);
end; % for iiW = 1:1:numel(www);

if SWT == 1;
figure('Name','Total Value in Conversion Region');
meshc(RRR, WWW, FF1); grid on; 
ylabel('Technology Level, W');  xlabel('Reserved Environment, R');
title(['Total Value in Conversion Region: \eta=', num2str(eta)], ...
       'FontSize', 14); 
else % if SWT == 1;
    tmpS=sprintf('Fig5Con_A%02.0fS%02.0fE%02.0f', tmpN);
    save(tmpS, 'RRR',      'WWW',      'FF1',      'KeyParam');
end; % if SWT == 1;

if SWT == 1;
figure('Name','Conversion Cost in Conversion Region');
meshc(RRR, WWW, cc); grid on; 
ylabel('Technology Level, W');  xlabel('Reserved Environment, R');
title(['Conversion Cost in Conversion Region: \eta=', num2str(eta)], ...
       'FontSize', 14); 
else % if SWT == 1;
    tmpS=sprintf('Fig6Cst_A%02.0fS%02.0fE%02.0f', tmpN);
    save(tmpS, 'RRR',      'WWW',      'cc',      'KeyParam');
end; % if SWT == 1;

%[C6] Total value in Combined Regions;
FFF     = FF1;                                % default = FF1;
for iiW = 1:1:numel(www);
    for iiR = 1:1:numel(rrr);
        if rrr(iiR) <= ppval(www(iiW), barr); % if R < R*, then
            FFF(iiW,iiR)   = FF0(iiW,iiR);    % FFF = FF0, otherwise FF1;
        end; % if rrr(iiR) <= ppval(www(iiW), Barr);
    end; % for iiR = 1:1:numel(rrr);
end; % for iiW = 1:1:numel(www);

if SWT == 1;
figure('Name','Total Value in Combined Regions');
meshc(RRR, WWW, FFF); grid on; 
ylabel('Technology Level, W');  xlabel('Reserved Environment, R');
title(['Total Value in Combined Regions: \eta=', num2str(eta)], ...
       'FontSize', 14); 
else % if SWT == 1;
    tmpS=sprintf('Fig7Tot_A%02.0fS%02.0fE%02.0f', tmpN);
    save(tmpS, 'RRR',      'WWW',      'FFF',      'KeyParam');
end; % if SWT == 1;


%% D: Equivalent Consumption Measure;
%[D1] Value alog Optimal Barrier;
%   => Recycle f0_ in [C5].     NB: f0_ is a function of W;

%[D2] Value alog Myopic Barrier 
%   => parallel to f0_ in [C5]. NB: again as a function of W;
if 324 == 324;                          % To be Not Used...
    br_my   = 0.0*ones(size(rr2));
    warning('off', 'MATLAB:quad:MaxFcnCount'     );
    warning('off', 'MATLAB:quad:ImproperFcnValue'); % disp(num2str(iiR));
for iiW = 1:1:numel(rr2);               % loop for integration
    br_my(iiW)= quad(Bprim, 0, rr2(iiW));% quad: adaptive simpson's rule;
end; % for iiW = 2:1:numel(rr2);
    warning('on',  'MATLAB:quad:MaxFcnCount'     );
    warning('on',  'MATLAB:quad:ImproperFcnValue');
br_my   = ph3*br_my;                    % B(R) along myopic barrier;
opVa_my = (  zzz.^bet  ).*( br_my );    % 3rd term along myopic barrier;
opVa_dy = (  zzz.^bet  ).*( bbr   );    % 3rd term along dynamic barrier;

f0_my   = zzz.*(1-rr2).^(1-1/eta)/(  (rho-alp)*(1-1/eta)  ) ...
        + rr2.^(1-1/eta)*phi/(   rho*(1-1/eta)  )           ...
        + opVa_my;                      % total value along the barrier;
end; % if 324 == 324;                   % To be Not Used...

%[D3] Equivalent Consumption Measure;
tmpCE       = NaN*ones(4,numel(www));
tmpCE(:,1)  =[0.0; NaN; NaN; NaN];
    warning('off', 'MATLAB:quad:MaxFcnCount'     );
    warning('off', 'MATLAB:quad:ImproperFcnValue');
for iiW = 2:1:numel(rr2);               % loop for each (Z, R(Z));
f0_CE   =@(Cadd_) ...
          zzz(iiW).*(1-rr2(iiW)+Cadd_).^(1-1/eta)/(rho-alp)/(1-1/eta)  ...
        + rr2(iiW).^(1-1/eta)*phi/(   rho*(1-1/eta)  )                 ...
        + ph3.*(  zzz(iiW).^bet  ).*quad(Bprim, 0, rr2(iiW), 1e-12) ...
        - f0_(iiW);
Cadd    = fzero(  f0_CE  ,  (rr1(iiW)-rr2(iiW))  );
tmpCE(:,iiW) = [www(iiW); f0_CE(0.0); Cadd/www(iiW); Cadd/(1-rr2(iiW))*100];
end; % for iiW = 2:1:numel(rr2);
    warning('on',  'MATLAB:quad:MaxFcnCount'     );
    warning('on',  'MATLAB:quad:ImproperFcnValue');
tmpCE4  = tmpCE(4,:);
tmpML   = (rr1-rr2)./rr1*100;

if 361 == 360;
figure('Name','Equivalent Consumption Measure for Myopic Loss');
if 400 == 40;
subplot(2,2,1);
plotyy(www, tmpCE(2,:), www, tmpCE(3,:)); grid on;
xlabel('Technology Level, W');
legend('CE to Myopic Loss', 'Loss in VF (LHS) ');
% legend for plotyy is weird!!!; 'Loss in VF (LHS)' is for tmpCE(2,:);
subplot(2,2,2);
[tmpAX,tmpH1,tmpH2] = plotyy(rr1, tmpCE(2,:), rr1, tmpCE(3,:)); grid on;
xlabel('Reserved Environment, R');
legend('CE to Myopic Loss', 'Loss in VF (LHS) '); 
set(tmpAX(1), 'XLim', [0,1]); set(tmpAX(2), 'XLim', [0,1]);
end; % if 400 == 40;
subplot(2,2,1);
plot(www, (rr1-rr2)./rr1*100 ); grid on; 
xlabel('Technology Level, W');
title('Horizontal Gap btw Myopic R and Optimal R (as % of Optimal R)');
subplot(2,2,2);
plot(rr1, (rr1-rr2)./rr1*100 ); grid on; 
xlabel('Reserved Environment, R'); set(gca, 'XLim', [0,1]);
title('Horizontal Gap btw Myopic R and Optimal R (as % of Optimal R)');
subplot(2,2,3);
plot(www, tmpCE( 4 ,:)); grid on; 
xlabel('Technology Level, W');
title('CE to Myopic Loss (as % of Actual C)');
subplot(2,2,4);
plot(rr1, tmpCE( 4 ,:)); grid on; 
xlabel('Reserved Environment, R'); set(gca, 'XLim', [0,1]);
title('CE to Myopic Loss (as % of Actual C)'); 
end % if 361 == 361;

if SWT == 1;
figure('Name','Equivalent Consumption Measure for Myopic Loss');
[tmp tmp1 tmp2] = ...
plotyy(www, tmpCE( 4 ,:),  www, (rr1-rr2)./rr1*100); grid on; 
xlabel('Technology Level, W');
set(tmp1, 'LineWidth', 2); set(tmp2, 'LineStyle', '-.');
legend('Consumption Equivalence (as % of Actual C)', ...
       'Horizontal Gap btw R^m and R^* (as % of Optimal R)');
title(['Consumption Equivalence for Myopic Loss: \eta=', num2str(eta)], ...
       'FontSize', 14);
title('It seems that legends order depends on the Matlab versions?');
else % if SWT == 1;
    tmpS=sprintf('Fig8ECL_A%02.0fS%02.0fE%02.0f', tmpN);
    save(tmpS, 'www',      'tmpCE4',     'tmpML',    'KeyParam');
end; % if SWT == 1;





%% E: Verification;
%debug;


toc;

