clc,clear;
addpath('/mq/home/m1mpg00/manuel_projects/durbinkoopman_mga')
addpath('/mq/home/m1mpg00/manuel_projects/PanelGDPTrendCycle/Data/Cleaned')
datady = xlsread('dy2018Q2.xlsx','B2:BB52'); % 2005Q2 to 2018Q2 51xT
datadu = xlsread('dURate_Quarter.xlsx','B2:FN52'); % 1976Q2 to 2018Q2
WitY2 = xlsread('weights2018Q2.xlsx','B2:BC52'); % 2005Q1 to 2004Q4
WitY = WitY2(:,2:end); % 2005Q2 to 2018Q2 
T = size(datady,2);
Tu = size(datadu,2);
n = size(datady,1);

datadu = datadu(:,Tu-T+1:end);
data = [datady;datadu];

statess = {'Alabama','Alaska','Arizona','Arkansas','California','Colorado','Connecticut','Delaware',...
           'District of Columbia','Florida','Georgia','Hawaii','Idaho','Illinois','Indiana','Iowa',...
           'Kansas','Kentucky','Louisiana','Maine','Maryland','Massachusetts','Michigan','Minnesota',...
           'Mississippi','Missouri','Montana','Nebraska','Nevada','New Hampshire','New Jersey',...
           'New Mexico','New York','North Carolina','North Dakota','Ohio','Oklahoma','Oregon','Pennsylvania',...
           'Rhode Island','South Carolina','South Dakota','Tennessee','Texas','Utah','Vermont',...
           'Virginia','Washington','West Virginia','Wisconsin','Wyoming'}';
%% Initial parameters

param = [1.5;-0.6;
        0.8*ones(n,1);ones(n,1);
        ones(n,1);ones(n,1);
        -0.25*ones(n,1);-0.25*ones(n,1);
        zeros(n,1);ones(n,1);
        ones(n,1);ones(n,1)];

phi1 = param(1);
phi2 = param(2);
mui = param(3:n+3-1);
alphai = param(n+3:2*n+3-1);
deltayi = param(2*n+3:3*n+3-1);
deltaui = param(3*n+3:4*n+3-1);
theta1i = param(4*n+3:5*n+3-1);
theta2i = param(5*n+3:6*n+3-1);
rhoi = param(6*n+3:7*n+3-1);
sigmazetai2 = param(7*n+3:8*n+3-1);
sigmaetayi2 = param(8*n+3:9*n+3-1);
sigmaetaui2 = param(9*n+3:10*n+3-1);
sigmae2 = 1;
sigmaetay2 = 1;
sigmaetau2 = 1;

%% Priors
m_phi = [1.5;-0.6];
invV_phi = diag([1/0.5^2;1/0.2^2]);
m_rho = 0;
invV_rho = 1;
m_beta_y = [0.8;1;1];
invV_beta_y = diag([1/0.8^2;1/1^2;1/1^2]);
m_beta_u = [1;-0.4;-0.1];
invV_beta_u = diag([1/1^2;1/0.1^2;1/0.1^2]);

% a = 2; b = 0.25; % inverse gamma with mean 0.25 and infinite variance (80% common trend - 20% idiosyncratic)
% a_upsilon = 2; b_upsilon = 3.23; % inverse gamma with mean 3.23 and infinite variance (80% common cycle - 20% idiosyncratic)
a = 2; b = 1; 
a_upsilon = 2; b_upsilon = 1; 
% a = 2; b = 1; % inverse gamma with mean 1 and infinite variance
% a = 3; b = 2; % inverse gamma with mean and variance 1
% a = 2.0001; b = 0.010001; % inverse gamma with mean 0.01 and variance 1
% a = 2; b = 0.01; % inverse gamma with mean 0.01 and infinite variance
%%
tic
m = 100000;
paramdraw = [param zeros(length(param),m)];
cycles = zeros(T,n+1,m);
for j = 2:m+1
    param = paramdraw(:,j-1);
    %% Obtaining draws of the states using Durbin-Koopman (2002)
    xt = DKdrawsGDPUnempStateAR1(data,param,T,n);

    %% Obtaining draws of parameters of common cycle
    cyc = xt(1:3,:)';
    c_t = cyc(:,1);
    c_tm1 = cyc(:,2);
    c_tm2 = cyc(:,3);
    Yc = c_t;
    Xc = [c_tm1 c_tm2];
    phibar = (invV_phi + Xc'*Xc/sigmae2)\(invV_phi*m_phi + Xc'*Yc/sigmae2);
    Vphi = inv(invV_phi + Xc'*Xc/sigmae2);
    phi = rmvnrnd(phibar',Vphi,1,[1 1; -1 1; 0 -1; 0 1],[1;1;1;1]);
    phi1 = phi(1);
    phi2 = phi(2);
    
    %% Obtaining draws of state-level parameters

    % The order of the latent variables is:
    % c_t, c_{t-1}, c_{t-3}, eta_{yt}, eta_{ut}, nu_{it}, nu_{it-2}, nu_{it-3}
    k1 = 5; % There are 5 common components
    ks1 = 0;
    ks2 = 3;
    ks = ks1 + ks2; % There are 3 latent variables nu_it (t,t-1,t-2)
    k2 = ks*n;
    k = k1+k2; % There are 5 common factors and 3 latent variables for each of the 51 states

    etay = xt(4,:)';
    etau = xt(5,:)';
    dc_t = c_t - c_tm1;
    dc_tm1 = c_tm1 - c_tm2;

    for i = 1:n
        % Obtaining rho_i and sigmazeta2_i
        upsilon_i = xt(k1+ks1+1+ks*(i-1):k1+ks*i,:)';
        upsilon_it = upsilon_i(:,1);
        upsilon_itm1 = upsilon_i(:,2);
        upsilon_itm2 = upsilon_i(:,3);
        Yrho = upsilon_it;
        Xrho = upsilon_itm1;
        rhobar = (invV_rho + Xrho'*Xrho/sigmazetai2(i))\(invV_rho*m_rho + Xrho'*Yrho/sigmazetai2(i));
        Vrho = inv(invV_rho + Xrho'*Xrho/sigmazetai2(i));
        rho_s = trandn((-1-rhobar)/sqrt(Vrho),(1-rhobar)/sqrt(Vrho));
        rhoi(i) = rhobar + sqrt(Vrho)*rho_s;
        erhoi = Yrho-rhoi(i)*Xrho;
        sigmazetai2(i) = 1/gamrnd(a_upsilon+0.5*T,1/(b_upsilon+0.5*(erhoi'*erhoi)));

        % Obtaining parameters of the observation equations
        dyi  = datady(i,:)';
        dui  = datadu(i,:)';
        dupsilon_it = upsilon_it - upsilon_itm1;
        dupsilon_itm1 = upsilon_itm1 - upsilon_itm2;
        ci_t = alphai(i)*c_t + upsilon_it;
        dci_t = alphai(i)*dc_t + dupsilon_it;
        dci_tm1 = alphai(i)*dc_tm1 + dupsilon_itm1;
        
        % deltau, theta1i, theta2i, and sigmaetaui2
        y = dui;
        X = [etau dci_t dci_tm1];
        betabar = (invV_beta_u + X'*X/sigmaetaui2(i))\(invV_beta_u*m_beta_u + X'*y/sigmaetaui2(i));
        Vbetabar = inv(invV_beta_u + X'*X/sigmaetaui2(i));
        betadraw = rmvnrnd(betabar',Vbetabar,1,[-1 0 0],0)';
        deltaui(i) = betadraw(1);
        theta1i(i) = betadraw(2);
        theta2i(i) = betadraw(3);
        etaui = y - X*betadraw;
        sigmaetaui2(i) = 1/gamrnd(a+0.5*T,1/(b+0.5*(etaui'*etaui)));
        
        % mui, alphai, deltayi, and sigmaetayi2
        dytilde = dyi - dupsilon_it;
        dutilde = dui - deltaui(i)*etau - theta1i(i)*dupsilon_it - theta2i(i)*dupsilon_itm1;
        y = [dytilde;dutilde];
        X = [ones(T,1) dc_t etay;zeros(T,1) theta1i(i)*dc_t+theta2i(i)*dc_tm1 zeros(T,1)];
        Sigmainv = diag([ones(T,1)/sigmaetayi2(i);ones(T,1)/sigmaetaui2(i)]);
        betabar = (invV_beta_y + X'*Sigmainv*X)\(invV_beta_y*m_beta_y + X'*Sigmainv*y);
        Vbetabar = inv(invV_beta_y + X'*Sigmainv*X);
        if i == 5 % California
            betadraw = rmvnrnd(betabar',Vbetabar,1,[0 -1 0; 0 0 -1; -1 0 0],[0;0;0])';
        else
            betadraw = rmvnrnd(betabar',Vbetabar,1,[0 0 -1; -1 0 0],[0;0])';
        end
        mui(i) = betadraw(1);
        alphai(i) = betadraw(2);
        deltayi(i) = betadraw(3);
        etayi = dytilde - [ones(T,1) dc_t etay]*betadraw;
        sigmaetayi2(i) = 1/gamrnd(a+0.5*T,1/(b+0.5*(etayi'*etayi)));
    end

%%
    param = [phi1;phi2;
             mui;alphai;
             deltayi;deltaui;
             theta1i;theta2i;
             rhoi;sigmazetai2;
             sigmaetayi2;sigmaetaui2];
     
%%
    paramdraw(:,j) = param;
    alpha = WitY'*alphai;
    upsilon_t_all_i = xt(k1+ks1+1:ks:end,:)';
    upsilon_t = sum(WitY'.*upsilon_t_all_i,2);    
    cycles(:,:,j) = [alpha.*c_t+upsilon_t c_t*alphai'+upsilon_t_all_i];

    j
end
toc;

%%
paramGDPUnempSt = paramdraw(:,0.5*m+1:m);
cyclesGDPUnempSt = cycles(:,:,0.5*m+1:m);

 
mean(paramGDPUnempSt(1:2,:),2)
[cellstr(statess) num2cell(mean(paramGDPUnempSt(3:n+3-1,:),2)) num2cell(mean(paramGDPUnempSt(n+3:2*n+3-1,:),2)) ...
 num2cell(mean(paramGDPUnempSt(2*n+3:3*n+3-1,:),2)) num2cell(mean(paramGDPUnempSt(3*n+3:4*n+3-1,:),2)) ...
 num2cell(mean(paramGDPUnempSt(4*n+3:5*n+3-1,:),2)) num2cell(mean(paramGDPUnempSt(5*n+3:6*n+3-1,:),2))]
[cellstr(statess) num2cell(mean(paramGDPUnempSt(6*n+3:7*n+3-1,:),2)) num2cell(sqrt(mean(paramGDPUnempSt(7*n+3:8*n+3-1,:),2)))]
[cellstr(statess) num2cell(sqrt(mean(paramGDPUnempSt(8*n+3:9*n+3-1,:),2))) num2cell(sqrt(mean(paramGDPUnempSt(9*n+3:10*n+3-1,:),2)))]

figure1 = figure;
tt = (2005.25:0.25:2017)';
plot(tt,mean(cyclesGDPUnempSt(:,1,:),3));
%% saving results
param_init = paramGDPUnempSt(:,end);
save('param_init','param_init');
