% UC_output_gap.m

clear; clc; 
nsim = 10000;
burnin = 1000;
data_raw = load('USGDP_2015Q4.csv');
data = 100*log(data_raw);
y = data;
T = length(y);

    % prior
a0 = [750;750]; B0 = 100*eye(2);
phi0 = [1.3 -.7]'; iVphi = speye(2);
nu_sigc2 = 3; S_sigc2 = 1*(nu_sigc2-1);
sigtau2_ub = .01;

    % initialize for storeage
store_theta = zeros(nsim,6); % [phi, sigc2, sigtau2, tau0]
store_tau = zeros(nsim,T); 
store_mu = zeros(nsim,T);    % annualized trend growth

    % initialize the Markov chain
phi = [1.34 -.7]';
tau0 = [y(1) y(1)]'; % [tau_{0}, tau_{-1}]
sigc2 = .5;
sigtau2 = .001;

    % construct a few things
H2 = speye(T) - 2*sparse(2:T,1:(T-1),ones(1,T-1),T,T) ...
    + sparse(3:T,1:(T-2),ones(1,T-2),T,T);
H2H2 = H2'*H2;
Hphi = speye(T) - phi(1)*sparse(2:T,1:(T-1),ones(1,T-1),T,T) + ...
    - phi(2)*sparse(3:T,1:(T-2),ones(1,T-2),T,T);
Xtau0 = [(2:T+1)' -(1:T)'];
n_grid = 500;
count_phi = 0; 

for isim = 1:nsim+burnin     
        % sample tau  
    alp_tau = H2\[2*tau0(1)-tau0(2);-tau0(1);sparse(T-2,1)];    
    Ktau = H2H2/sigtau2 + Hphi'*Hphi/sigc2;
    tau_hat = Ktau\(H2H2*alp_tau/sigtau2 + Hphi'*Hphi*y/sigc2);
    tau = tau_hat + chol(Ktau,'lower')'\randn(T,1);

        % sample phi
    c = y-tau;
    Xphi = [[0;c(1:T-1)] [0;0;c(1:T-2)]];    
    Kphi = iVphi + Xphi'*Xphi/sigc2;
    phi_hat = Kphi\(iVphi*phi0 + Xphi'*c/sigc2);
    phic = phi_hat + chol(Kphi,'lower')'\randn(2,1);
    if sum(phic) < .99 && phic(2) - phic(1) < .99 && phic(2) > -.99
        phi = phic;    
        Hphi = speye(T) - phi(1)*sparse(2:T,1:(T-1),ones(1,T-1),T,T)...
            - phi(2)*sparse(3:T,1:(T-2),ones(1,T-2),T,T);
        count_phi = count_phi + 1;    
    end
    
        % sample sigc2
    sigc2 = 1/gamrnd(nu_sigc2 + T/2,1/(S_sigc2 ...
        + (c-Xphi*phi)'*(c-Xphi*phi)/2));
    
        % sample sigtau2
    del_tau = [tau0(1);tau(1:T)] - [tau0(2);tau0(1);tau(1:T-1)]; 
    f_tau = @(x) -T/2*log(x) ...
        - sum((del_tau(2:T) - del_tau(1:T-1)).^2)./(2*x);       
    sigtau2_grid = linspace(rand/1000,sigtau2_ub-rand/1000,n_grid);
    lp_sigtau2 = f_tau(sigtau2_grid);
    p_sigtau2 = exp(lp_sigtau2-max(lp_sigtau2));
    p_sigtau2 = p_sigtau2/sum(p_sigtau2);
    cdf_sigtau2 = cumsum(p_sigtau2);
    sigtau2 = sigtau2_grid(find(rand<cdf_sigtau2, 1 ));    
    
        % sample tau0 
    Ktau0 = B0\speye(2) + Xtau0'*H2H2*Xtau0/sigtau2;
    tau0_hat = Ktau0\(B0\a0 + Xtau0'*H2H2*tau/sigtau2);        
    tau0 = tau0_hat + chol(Ktau0,'lower')'\randn(2,1);    
    
    if (mod(isim, 1000) ==0)
        disp([num2str(isim) ' loops... ']);
    end
    
    if isim > burnin
        i = isim-burnin;
        store_tau(i,:) = tau';
        store_theta(i,:) = [phi' sigc2 sigtau2 tau0'];
        store_mu(i,:) = 4*(tau-[tau0(1);tau(1:end-1)])';
    end    
end
       
tau_hat = mean(store_tau)';
theta_hat = mean(store_theta)'
theta_CI = quantile(store_theta,[.025 .975])
mu_hat = mean(store_mu)';

%% plot of graphs
tt = (1947:.25:2015.75)';
figure;
hold on
    plot(tt,y-tau_hat,'linewidth',1);
    plot(tt,zeros(T,1),'--k','linewidth',1);
hold off
title('Output gap estimates');
xlim([1947 2016]); box off;
set(gcf,'Position',[100 100 800 300]);

figure;
plot(tt,mu_hat); xlim([1947 2016]); ylim([1 4.5]); box off;
set(gcf,'Position',[100 100 800 300]);
title('Trend output growth estimates');
