% Bayes Analsyis of d_tfp .. graphs 
% 
clear all;
small = 1.0e-10;
big = 1.0e+6;
this_date = datestr(now,'yyyymmdd');

% -----  File Directories -- overhead, etc.
 outdir = 'out/';
 figdir = 'fig/';
 matdir = 'mat/';
 rtsdir = '../rts/';
 p_initial = path;
 vardef_dir = 'variable_definitions/';            % Directory for necessary matlab_procs and functions                         
 p1=path(vardef_dir,p_initial);
 utility_dir = 'm_utilities';
 path(utility_dir,p1);

% -- Load Data
load_data=0;         % Set = 1 first time to read data from Excel files (then they are saved as MATLAB file)
fhsw_bpea_data;      % fhsw_bpea_data reads in the full dataset .. all variables, etc. saved in datain.xx

% Get Supply Variables
supply_variables;


% Cyclcical Adjustment Parameters
du_nlead = 2;                   % number of leads
du_nlag = 2;                    % number of lags
  
% -- Set up regressors for cylical adjustment
d_urate = dif(cps_urate,1);     % first difference of unemployment rate
xca = lag(d_urate,-du_nlead);   % xca is cyclical adjustment variable
for i = (-du_nlead+1):du_nlag;
 xca = [xca lag(d_urate,i)];
end;

% ------------------ Analysis for 1956-2016 Sample Period --------------
rng(996530);
ismpl_ol = smpl(calvec,[1981 3],[2016 2],4);  % Sample period for estimating Cyclical adjustment coefficients
dy = d_bus_tfp;
bw_bw = 60;
[dy_cycle,dy_ca,dy_mu_60,dy_irregular,var_dy_mu_60,b_okun_sum,se_b_okun_sum ] = cyclical_adjustment(dy,xca,ismpl_ol,bw_bw);
dy4_ca = (dy_ca + lag(dy_ca,1) + lag(dy_ca,2) + lag(dy_ca,3))/4;
se_dy_mu_60 = sqrt(diag(var_dy_mu_60));


% RW + WN model for dy_ca;
% y(t) = tau(t) + eps(t); tau(t)=tau(t-1)+eta(t)
% Data
y = dy_ca(ismpl_ol==1);
tmp = (dy_ca + lag(dy_ca,1) + lag(dy_ca,2) + lag(dy_ca,3))/4;
y4 = tmp(ismpl_ol==1);
T = size(y,1);

% Prior for tau(0)
mu_0 = 1;
var_0 = 10;  % Keep This fixed;

% Prior on Standard Deviation of Eps -- uniform on sd_eps_min to sd_eps_max;
sd_y = std(y);
sd_eps_min = sd_y*(2/3);
sd_eps_max = sd_y*(4/3);
n_sd_eps = 20;
sd_eps_values = linspace(sd_eps_min,sd_eps_max,n_sd_eps)';

% prior on standard deviation of eta
sd_eta_min = 0.001;
sd_4y_change_max = 1;  % largest value for sd of chance in tau over 4-year period
%sd_4y_change_max = .001;
sd_eta_max = sd_4y_change_max/4;  % 4 = sqrt(16 quarters)
n_sd_eta = 40;
sd_eta_values = linspace(sd_eta_min,sd_eta_max,n_sd_eta)';

% RW Variance
tmp = tril(ones(T,T));
rw_var = tmp*tmp';

% Variance from initial values
var_0_mat = var_0*ones(T,T);
y_m0 = y - mu_0;

% Construct Log-likelihood values -- brute force
nvalues = n_sd_eps*n_sd_eta;
llf = NaN(nvalues,3);
n=0;
for i = 1:n_sd_eps;
    for j = 1:n_sd_eta;
        n=n+1;
        llf(n,1) = sd_eps_values(i);
        llf(n,2) = sd_eta_values(j);
    end;
end;

tau_mean_mat = NaN(T,nvalues);
tau_var_mat = NaN(T,T,nvalues);
cholesky_tau_var_mat = NaN(T,T,nvalues);
tausq_mean_mat = NaN(T,nvalues);
for n = 1:nvalues;
    var_eps = llf(n,1)^2;
    var_eta = llf(n,2)^2;
    v_tau = var_eta*rw_var + var_0_mat;
    v_eps = var_eps*eye(T);
    sig = v_tau+v_eps;
    det_sig = det(sig);
    sig_inv = inv(sig);
    qf = y_m0'*sig_inv*y_m0;
    llf(n,3) = -0.5*(log(det_sig)+qf);
    K = v_tau*sig_inv;
    tau_mean = mu_0 + K*(y_m0);
    tau_var = (eye(T)-K)*v_tau;
    tau_mean_mat(:,n) = tau_mean;   
    tau_var_mat(:,:,n)=tau_var;
    cholesky_tau_var_mat(:,:,n)=chol(tau_var)';
end;
lf = llf;
lf(:,3) = exp(lf(:,3));
post = lf(:,3)/sum(lf(:,3));
post = [lf(:,1:2) post];

% Posterior probabilities
post_eta = NaN(n_sd_eta);
prior_eta = ones(n_sd_eta)/n_sd_eta;
for i = 1:n_sd_eta;
    ii = post(:,2)==sd_eta_values(i);
    post_eta(i) = sum(post(ii==1,3));
end;

tau_mean = NaN(dnobs,1);
tau_mean(ismpl_ol == 1) = tau_mean_mat*post(:,3);

% Get draws of Y
ndraw = 10000;
postvec = cumsum(post(:,3));
calvec_ol = calvec(ismpl_ol==1);
ismpl = smpl(calvec_ol,[1981 3],[2016 2],4);
cal = calvec_ol(ismpl==1);
iT = find(cal==(2007 + (4-1)/4));
stemp = NaN(ndraw,4);
tau_draw = NaN(T,ndraw);
for i = 1:ndraw;
    u = rand(1);
    il = find(cumsum(postvec >= u) == 1);
    m = tau_mean_mat(:,il);
    c = squeeze(cholesky_tau_var_mat(:,:,il));
    ydraw = m + c*randn(size(m,1),1);
    tau_draw(:,i) = ydraw;
    yy = ydraw(ismpl==1);
    [maxy,ii]=max(yy);
    d1981 = maxy-yy(1);
    d2007 = maxy-yy(iT);
    d2007_2016 = yy(iT)-yy(end);
    stemp(i,:) = [cal(ii) d1981 d2007 d2007_2016];
end;

pct = [0.16666 0.50 0.83333]';

% Posterior probabilities
tau_pct = NaN(T,3);
for t = 1:T;
    tau_pct(t,:)=pctile(tau_draw(t,:)',pct);
end;


ismpl1 = smpl(calvec,[1981 3],[1995 4],4);
ismpl2 = smpl(calvec,[1996 1],[2006 1],4);
ismpl3 = smpl(calvec,[2006 2],[2016 2],4);
dy_mu_break = NaN(dnobs,1);
dy_mu_break(ismpl1==1) = mean(dy_ca(ismpl1==1));
dy_mu_break(ismpl2==1) = mean(dy_ca(ismpl2==1));
dy_mu_break(ismpl3==1) = mean(dy_ca(ismpl3==1));

figure;
axis_font_size = 20;
tmp = NaN(dnobs,3);
tmp(ismpl_ol==1,:)=tau_pct;
plot(calvec(ismpl_ol==1),dy4_ca(ismpl_ol==1),'- k','LineWidth',1);
hold on;
 plot(calvec(ismpl_ol==1),dy_mu_60(ismpl_ol==1),' - b','LineWidth',3);
 plot(calvec(ismpl_ol==1),dy_mu_break(ismpl_ol==1),' - k','LineWidth',3);
 cal = calvec(ismpl_ol==1);
 xu = tmp((ismpl_ol==1),3);
 xl = tmp((ismpl_ol==1),1);
 for t = 1:size(cal,1)-1;
  patch([cal(t) cal(t) cal(t+1) cal(t+1)],[xl(t) xu(t) xu(t) xl(t)],'k','EdgeColor','none');
  alpha(0.1);
 end;
hold off;
xlim([1980 2020]);
ax = gca;  % This gets axis limits for current plot (must be open)
ax.FontSize = axis_font_size;
savefig([figdir 'Fig_8_' this_date]);
saveas(gcf,[figdir 'Fig_8.epsc']);


path(p_initial);  % Reset path