%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%     TASK 4


%%%     c ~ N(nu, sigma)


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%





clear all;
close all;





%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%


%%%%%%   SET PARAMETER VALUES HERE


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%

s = 2;
nu = 30;
sigma = 7;

%%  This is a huge loop that varies s, nu and sigma for the whole code!


while s<=2
   
   while nu<=30
      
      while sigma <= 7
         
         
         
     
%%%%% VISUALIZE THE FUNCTION FOR AN INFORMED GUESS OF THE INITIAL VALUES
%c_try = nu-3*sigma:5:nu+3*sigma;
%u = size(c_try);
%u = u(1,2)
%y = [];
%for i=1:u
%y(i) = fun_exp(c_try(i), s, nu, sigma);
%i
%end
%plot(c_try,y)





%%%  First, compute C_star for exponential discounter


%%%  fun_exp is a function whose output is s - c_star + theta, where
%%%  s is the daily cost of default strategy, c is the switching cost, and theta is EV.
%%%  c* is the c value that makes fun_exp = 0




%if (fun_exp(nu-3*sigma, s, nu, sigma)>0) + (fun_exp(nu+3*sigma, s, nu, sigma)>0) == 1
   [c_star, fval] = fzero('fun_exp', [nu-3*sigma, nu+3*sigma], [], s, nu, sigma);
%else
 %  c_star = 0;
%end




%%% Try different values of beta.

beta = 0.5:0.1:1;
b = size(beta);
b = b(1,2);



%%% C_star for naive discounter is equal to beta * exponential C*
%%%  ED = Expected number of failures before success  formula = (1-F) / F
%%%%  P120 is the probability that there are at least 120 failures for the first 120 days = (1-F)^120

c_naive = beta  .* c_star;

F_naive = normcdf(c_naive, nu, sigma);

ED_naive = (1-F_naive) ./ F_naive;

P120_naive = (1-F_naive).^120;



%%% C*  Sophisticated:  Follow Laibson's appendix.

c_sop = [];


%%% for each beta, find C* sophisticated by finding the c value that makes fun_sop = 0
%%%  fun_sop is a function whose output is beta*(s + theta) - c* where
%%%  s is the daily cost of default strategy, c is the switching cost, and theta is EV.


for i=1:b
   
   multiplier = 3;
   %if sigma >=5 multiplier = 1.5; end
   %if fun_sop(nu-multiplier*sigma, s, nu, sigma, beta(i))>0 + fun_sop(nu+multiplier*sigma, s, nu, sigma,beta(i))>0 == 1
      c_sop(i) = fzero('fun_sop', [nu-multiplier*sigma, nu+multiplier*sigma], [], s, nu, sigma,beta(i));
   %else
   %  c_sop(i) = 0;
   %end
   

end


F_sop = normcdf(c_sop, nu, sigma);


ED_sop = (1-F_sop) ./ F_sop;


P120_sop = (1-F_sop).^120;




figure;

subplot(2,2,1), plot(beta,c_naive, '.', beta, c_sop, '-');
title(['C*, nu = ',num2str(nu),' sigma = ', num2str(sigma), ' s = ', num2str(s)]);
xlabel('beta');
ylabel('C*');
%axis([0.5 1 0 30]);
legend('C* naive', 'C* sop', 4);

subplot(2,2,2),plot(beta,F_naive, '.', beta,F_sop,'-');
title(['F(C*), nu = ',num2str(nu),' sigma = ', num2str(sigma), ' s = ', num2str(s)]);
xlabel('beta');
ylabel('F(C*)');
axis([0.5 1 0 1]);
legend('F(C*) naive', 'F(C*) sop', 4);

subplot(2,2,3),plot(beta,ED_naive, '.', beta, ED_sop,'-');
title(['E[days], nu = ',num2str(nu),' sigma = ', num2str(sigma), ' s = ', num2str(s)]);
xlabel('beta');
ylabel('ED');
axis([0.5 1 0 100]);
legend('ED naive', 'ED sop', 4);

subplot(2,2,4),plot(beta,P120_naive, '.', beta,P120_sop,'-');
title(['P120, nu = ',num2str(nu),' sigma = ', num2str(sigma), ' s = ', num2str(s)]);
xlabel('beta');
ylabel('P120');
axis([0.5 1 0 0.7]);
legend('P120 naive', 'P120 sop', 4);





      if sigma == 1
         sigma = 2;
      elseif sigma == 2
         sigma = 5;
      elseif sigma == 5
         sigma = 7;
      else
         sigma = 9;
      end
      
    %% this end is for the innermost loop for sigma  
    s
    nu
    sigma
   
 end
 
   sigma =1;  
   if nu == 10 
      nu = 15;
   elseif nu == 15
      nu = 30;
   else nu = 31;
      
   end
   
   
   s
   nu
   sigma

   %% this end is for the nu loop  nu ~ {10, 15, 30}
end

nu = 10;

s = s+1;

s
nu
sigma


%% this end is for the s loop, the outermost one.  s ~ {1,2}
end