% VAR using the Gibbs sampler, based on independent Normal Wishart prior

%--------------------------------------------------------------------------
clear all;
clc;
addpath('function')

%------------------------------LOAD DATA-----------------------------------
dataset = 'italy_euro_oo_sr'; % Type in the data set to be used

all_var= xlsread(dataset);
 t_code = xlsread(['t_code_' dataset]);
[T, K] = size(all_var);

% Transform data
YYraw = zeros(T,K);

for i=1:K
 YYraw(:,i) = transx(all_var(:,i),t_code(i,1));
end

YYraw=YYraw;

%----------------------------PRELIMINARIES---------------------------------
% Define specification of the VAR model
constant = 1;        % 1: if you desire intercepts, 0: otherwise 
p = 3;               % Number of lags on dependent variables                    
stability =1;	     % Checks whether draws are stable

% Set prior for BVAR model:
prior = 2;  % prior = 1 --> Indepependent Normal-Whishart Prior
            % prior = 2 --> Indepependent Minnesota-Whishart Prior

% Gibbs-related preliminaries
nsave = 1000;         % Final number of draws to save
nburn = 2000;         % Draws to discard (burn-in)
ntot = nsave + nburn;  % Total number of draws
it_print = 100;      % Print on the screen every "it_print"-th iteration
offset = 0.0000000001;

%------------------------------IRF-horizon---------------------------------
hor=40;

%------------- Help: determine number of variables and number of parameters
[T_help nvar_help] = size(YYraw);

if constant==1
    K_help=nvar_help*p+1;
else
    K_help=nvar_help*p;
end


%------------------------------Set array of window size--------------------
% Please note, that futher changes in the code are necessary to obtain the
% full sample estimate. You should not set linespace to the full sample!
window_size= linspace(96+p,96+p,1);

%determine number of rolling windows
n_roll = T_help-max(window_size)+1;

%------------------------------generate empty matrices---------------------


ALPHA_draws_win_IF = zeros(K_help, nvar_help,  n_roll, size(window_size,2));
SIGMA_draws_win_IF = zeros(nvar_help, nvar_help, n_roll, size(window_size,2));

% matrices to store irf
% irf_all_win = zeros(nvar_help,nvar_help,hor,nsave,n_roll,size(window_size,2));
irf_median_win = zeros(nvar_help,nvar_help,hor, n_roll,size(window_size,2));
irf_lb_win = zeros(nvar_help,nvar_help,hor, n_roll,size(window_size,2));
irf_ub_win = zeros(nvar_help,nvar_help,hor, n_roll,size(window_size,2));

% save alpha and sigma
ALPHA_draws_win=zeros(K_help,nvar_help,1000,n_roll,size(window_size,2));
SIGMA_draws_win=zeros(nvar_help,nvar_help,1000,n_roll,size(window_size,2));

tic
for iii =1:size(window_size,2) % increases the window size (72,73,74,...,120)
    
    for ii=1:n_roll % moves the window forward (adds another observation)
        if n_roll ==1
            Yraw=YYraw;
        else
            Yraw = YYraw((max(window_size)-min(window_size)-iii+ii+1):(max(window_size)-1+ii),:);
        end
        
%--------------------------DATA HANDLING-----------------------------------
% Get initial dimensions of dependent variable
[Traw M] = size(Yraw);

% Rename variables (due to manipulation of original code)
Y1 = Yraw;
Y2 = Yraw;

        
% Generate lagged Y matrix. This will be part of the X matrix
Ylag = mlag2(Y2,p); % Y is [T x M]. ylag is [T x (Mp)]

% Now define matrix X which has all the R.H.S. variables (constant, lags of
% the dependent variable and exogenous regressors/dummies)
if constant==1
    X1 = [ones(Traw-p,1) Ylag(p+1:Traw,:)];
else
    X1 = Ylag(p+1:Traw,:);  %#ok<UNRCH>
end

% Get size of final matrix X
[Traw3 K] = size(X1);

% Create the block diagonal matrix Z
Z1 = kron(eye(M),X1);

% Form Y matrix accordingly
% Delete first "LAGS" rows to match the dimensions of X matrix
Y1 = Y1(p+1:Traw,:); 

% Traw was the dimesnion of the initial data. T is the number of actual 
% time series observations of Y and X (we lose the p-lags)
T = Traw - p;

% Rename variables (due to manipulation of original code)
Y = Y1;
X = X1;
Z = Z1;



%-----------------------------PRELIMINARIES--------------------------------
% First get ML estimators
A_OLS = inv(X'*X)*(X'*Y); % This is the matrix of regression coefficients
a_OLS = A_OLS(:);         % This is the vector of parameters, i.e. it holds
                          % that a_OLS = vec(A_OLS)
SSE = (Y - X*A_OLS)'*(Y - X*A_OLS);   % Sum of squared errors
SIGMA_OLS = SSE./(T-K+1);

% Initialize Bayesian estimation
alpha = a_OLS;     % This is the single draw from the posterior of alpha
ALPHA = A_OLS;     % This is the single draw from the posterior of ALPHA
SSE_Gibbs = SSE;   % This is the single draw from the posterior of SSE
SIGMA = SIGMA_OLS; % This is the single draw from the posterior of SIGMA

% Storage space for posterior draws
alpha_draws = zeros(K*M,nsave);
ALPHA_draws = zeros(K,M,nsave);
SIGMA_draws = zeros(M,M,nsave);

%-----------------Prior hyperparameters for bvar model
n = K*M; % Total number of parameters (size of vector alpha)
% Define hyperparameters
if prior == 1 % Normal-Wishart
    A_prior = [zeros(1,M); 1*eye(M); zeros((p-1)*M,M)];  %<---- prior mean of ALPHA (parameter matrix) 
    a_prior = A_prior(:); 
    V_prior = 10*eye(n);     %<---- prior variance of alpha
    
    % Hyperparameters on inv(SIGMA) ~ W(v_prior,inv(S_prior))
    v_prior = M+1;             %<---- prior Degrees of Freedom (DoF) of SIGMA
    S_prior = eye(M);            %<---- prior scale of SIGMA
    inv_S_prior = inv(S_prior);
elseif prior == 2 % Minnesota-Whishart
    % Prior mean on VAR regression coefficients
    A_prior = [zeros(1,M); 1*eye(M); zeros((p-1)*M,M)];  %<---- prior mean of ALPHA (parameter matrix) 
    a_prior = A_prior(:);               %<---- prior mean of alpha (parameter vector)
    
    % Minnesota Variance on VAR regression coefficients
    % First define the hyperparameters 'a_bar_i'
    a_bar_1 = 0.5;
    a_bar_2 = 0.25;
    a_bar_3 = 10^2;
    
    % Now get residual variances of univariate p-lag autoregressions. Here
    % we just run the AR(p) model on each equation, ignoring the constant
    % and exogenous variables (if they have been specified for the original
    % VAR model)
    sigma_sq = zeros(M,1); % vector to store residual variances
    for i = 1:M
        % Create lags of dependent variable in i-th equation
        Ylag_i = mlag2(Yraw(:,i),p);
        Ylag_i = Ylag_i(p+1:Traw,:);
        % Dependent variable in i-th equation
        Y_i = Yraw(p+1:Traw,i);
        % OLS estimates of i-th equation
        alpha_i = inv(Ylag_i'*Ylag_i)*(Ylag_i'*Y_i);
        sigma_sq(i,1) = (1./(Traw-p+1))*(Y_i - Ylag_i*alpha_i)'*(Y_i - Ylag_i*alpha_i);
    end
    
    % Now define prior hyperparameters.
    % Create an array of dimensions K x M, which will contain the K diagonal
    % elements of the covariance matrix, in each of the M equations.
    V_i = zeros(K,M);
    
    % index in each equation which are the own lags
    ind = zeros(M,p);
    for i=1:M
        ind(i,:) = constant+i:M:K;
    end
    for i = 1:M  % for each i-th equation
        for j = 1:K   % for each j-th RHS variable
            if constant==1 % if there is constant in the model use this code:
                if j==1
                    V_i(j,i) = a_bar_3*sigma_sq(i,1); % variance on constant                
                elseif find(j==ind(i,:))>0
                    V_i(j,i) = a_bar_1./(p^2); % variance on own lags           
                else
                    for kj=1:M
                        if find(j==ind(kj,:))>0
                            ll = kj;                   
                        end
                    end                 % variance on other lags   
                    V_i(j,i) = (a_bar_2*sigma_sq(i,1))./((p^2)*sigma_sq(ll,1));      
                end
            else  % if there is no constant in the model use this:
                if find(j==ind(i,:))>0
                    V_i(j,i) = a_bar_1./(p^2); % variance on own lags
                else
                    for kj=1:M
                        if find(j==ind(kj,:))>0
                            ll = kj;
                        end                        
                    end                 % variance on other lags  
                    V_i(j,i) = (a_bar_2*sigma_sq(i,1))./((p^2)*sigma_sq(ll,1));            
                end
            end
        end
    end
    
    % Now V is a diagonal matrix with diagonal elements the V_i
    V_prior = diag(V_i(:));  % this is the prior variance of the vector alpha
    
    % Hyperparameters on inv(SIGMA) ~ W(v_prior,inv(S_prior))
    v_prior = M+1;
    S_prior = diag(sigma_sq); %eye(M);
    inv_S_prior = inv(S_prior);   
end
    
%========================== Start Sampling ================================
%==========================================================================
%tic;
%disp('Number of iterations');
for irep = 1:ntot  %Start the Gibbs "loop"
    if mod(irep,it_print) == 0
%         disp(irep);
%         toc;
    end
    
    VARIANCE = kron(inv(SIGMA),speye(T));
    V_post = inv(inv(V_prior) + Z'*VARIANCE*Z); % see Koop and Korobilis: Bayesian Multivariate Time Series for empirical Macroeconomic, S.13
    a_post = V_post*(inv(V_prior)*a_prior + Z'*VARIANCE*Y(:));
    
    % draw autoregressive coefficients and check for stability
    if stability ==1
        check = -1;
        while check<0
            try
            alpha = a_post + chol(V_post+offset)'*randn(n,1); % Draw of alpha
            catch
            alpha =alpha;   
            end
            ALPHA = reshape(alpha,K,M); % Draw of ALPHA
            stab = stability_BVAR(ALPHA,p,M);

            if stab==0
                check=10;
            end
        end
    else
        alpha = a_post + chol(V_post+offset)'*randn(n,1); % Draw of alpha
        ALPHA = reshape(alpha,K,M); % Draw of ALPHA
    end
    
    % Posterior of SIGMA|ALPHA,Data ~ iW(inv(S_post),v_post)
    v_post = T + v_prior;
    S_post = S_prior + (Y - X*ALPHA)'*(Y - X*ALPHA);
    
    SIGMA = inv(wish(inv(S_post),v_post));% Draw SIGMA    
    
    % Store results  
    if irep > nburn               
               
        %----- Save draws of the parameters
        alpha_draws(:,irep-nburn) = alpha;
        ALPHA_draws(:,:,irep-nburn) = ALPHA;
        SIGMA_draws(:,:,irep-nburn) = SIGMA;

    end % end saving results
       
end %end the main Gibbs for loop

%====================== End Sampling Posteriors ===========================
%toc;
% NOTE: Due to memory considerations I can't save the draws for alpha and sigma!
for j=1:K
    for jj=1:M
        test= momentg(squeeze(ALPHA_draws(j,jj,:)));        
        ALPHA_draws_win_IF(j,jj,ii,iii)=1/test.rne1; % 1/ because IF is inverse of rne, see primiceri 2005 Appendix
    end
end

for j=1:M
    for jj=1:M
        test= momentg(squeeze(SIGMA_draws(j,jj,:)));        
        SIGMA_draws_win_IF(j,jj,ii,iii)=1/test.rne1; % 1/ because IF is inverse of rne, see primiceri 2005 Appendix
    end
end

ALPHA_draws_win(:,:,:,ii,iii)=ALPHA_draws;
SIGMA_draws_win(:,:,:,ii,iii)=SIGMA_draws;

[irf_mat_all, irf_median, irf_ub, irf_lb,fevd_median, fevd_ub, fevd_lb, total_spillover_median, total_spillover_ub, total_spillover_lb, spillover_from_median, spillover_from_ub, spillover_from_lb, spillover_to_median, spillover_to_lb, spillover_to_ub] = IRF1(SIGMA_draws, ALPHA_draws, hor,  0.95, 0.05, 0, 0, 0);

% matrices to store irf
%irf_all_win(:,:,:,:,ii,iii)=irf_mat_all;
irf_median_win(:,:,:,ii,iii)=irf_median;
irf_lb_win(:,:,:,ii,iii)=irf_lb;
irf_ub_win(:,:,:,ii,iii)=irf_ub;

    disp(['Running the ', num2str(ii), 'of', num2str(n_roll), 'Rolling Windows'] )
    end
    disp(iii)
    toc
end
toc

name = [dataset, ', lag=' num2str(p) ', prior=' num2str(prior) ', stability=' num2str(stability) ', windowsize=96, ki=90, priormean=1'];
save(name) 
clear