close;
clear; 
clc;

%Step 1: Set parameter according to p 877
durgd=8;  %duration of good spell
durbd=8;  %duration of bad spell
durug=1.5; %duration of unemp in good state
durub=2.5; %duration of unemp in bad state
unempg=.04;  %unemployment rate in good times
unempb=.1;   %unemployment rate in bad times

      pgg00 = (durug-1)/durug;  %defined by parameters
      pbb00 = (durub-1)/durub;  %defined by parameters
      pbg00 = 1.25*pbb00;  %restriction from p 877
      pgb00 = 0.75*pgg00;  %restriction from p 877
      pgg01 = (unempg - unempg*pgg00)/(1-unempg); %restriction p872
      pbb01 = (unempb - unempb*pbb00)/(1-unempb);
      pbg01 = (unempb - unempg*pbg00)/(1-unempg);
      pgb01 = (unempg - unempb*pgb00)/(1-unempb);
      pgg = (durgd-1)/durgd;   %defined by parameters
      pgb = 1 - (durbd-1)/durbd;

      pgg10 = 1 - (durug-1)/durug; %defined by parameters
      pbb10 = 1 - (durub-1)/durub; %defined by parameters
      pbg10 = 1 - 1.25*pbb00;
      pgb10 = 1 - 0.75*pgg00;
      pgg11 = 1 - (unempg - unempg*pgg00)/(1-unempg); %restriction p872
      pbb11 = 1 - (unempb - unempb*pbb00)/(1-unempb);
      pbg11 = 1 - (unempb - unempg*pbg00)/(1-unempg);
      pgb11 = 1 - (unempg - unempb*pgb00)/(1-unempb);
      pbg = 1 - (durgd-1)/durgd;
      pbb = (durbd-1)/durbd;
pr=zeros(4,4);
      pr(1,1) = pgg*pgg11;  
      pr(2,1) = pbg*pbg11;
      pr(3,1) = pgg*pgg01;
      pr(4,1) = pbg*pbg01;
      pr(1,2) = pgb*pgb11;
      pr(2,2) = pbb*pbb11;
      pr(3,2) = pgb*pgb01;
      pr(4,2) = pbb*pbb01;
      pr(1,3) = pgg*pgg10;
      pr(2,3) = pbg*pbg10;
      pr(3,3) = pgg*pgg00;
      pr(4,3) = pbg*pbg00;
      pr(1,4) = pgb*pgb10;
      pr(2,4) = pbb*pbb10;
      pr(3,4) = pgb*pgb00;
      pr(4,4) = pbb*pbb00;
      
pr;
pr_sta=pr^200;
%calculate the eigenvalues and vectors of the stationary trans matrix
[v,d]=eig(pr_sta);
%this is the eigenvector that corresponds to the unit eigenvalue. 
pr_sta_unitevec = v(:,2)
%check that this is indeed the correct eigenvector. check_vec should be
%equal to zero.
check_evec = (pr_sta-eye(4,4))*pr_sta_unitevec;
%now normalize the sum of the elements of the eigen vector to one
pr_sta_unitevec_norm = pr_sta_unitevec./sum(pr_sta_unitevec,1);
disp(['Joint Asymp. Stationary Distribution = [',num2str(pr_sta_unitevec_norm'),']'])

%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%This is the second part of Rudi's handout for point, find the agg.
%transition prob.
pr_agg_pose = [pgg pgb; pbg pbb];
pr_agg_sta = (pr_agg_pose')^200;
[v,d]=eig(pr_agg_sta);
%this is the eigenvector that corresponds to the unit eigenvalue. 
pr_agg_sta_unitevec = v(:,2)
%check that this is indeed the correct eigenvector. check_vec should be
%equal to zero.
check_evec = (pr_agg_sta-eye(2,2))*pr_agg_sta_unitevec
%now normalize the sum of the elements of the eigenvector to one
pr_agg_sta_unitevec_norm = pr_agg_sta_unitevec./sum(pr_agg_sta_unitevec,1);
disp(['Agg. Asymp. Stationary Distribution = [',num2str(pr_agg_sta_unitevec_norm'),']'])


%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
%See the appendix of KS p. 891 for definition of the (k,kbar) space.
%Make kbar equidistant and fairly far spaced
kbar = [0:0.5:30]';
%Make k not equidistant and lots of observations near zero. The quadratic
%relationship just does this quickly.
k =  ([0:0.1:3]').^3;

% 
% 
% rand('seed',123456);
% share_g = rand(1,1);
% 
% %Put share_g in first column
% T = 1200;
% agg_share = 99*ones(2,T);
% agg_share(:,1) = [share_g; 1-share_g];
% for t = 2:1:T
%     agg_share(:,t) = pr_agg_pose'*agg_share(:,t-1);
% end
% 
% pr_agg_sta = (pr_agg_pose')^300;
% 
% for t = 2:1:T
%     agg_share(:,t) = pr_agg_sta*agg_share(:,t-1);
% end
% 
% % temp = rand(1200,1);
% % temp_array = temp
% % agg_1200 = pr_agg_sta*
% % agg_1200 = 0.5*temp+0.5*(1-temp);
% 
%