function [A,SIGMA,U,V] = bvarc(Yraw,p,prior)

% Y are the VAR data, p is number of lags
% prior = 1 --> Noninformative Prior
% prior = 2 --> Minnesota Prior
% prior = 3 --> Natural conjugate Prior
  
%--------------------------DATA HANDLING-----------------------------------
% Get initial dimensions of dependent variable
[Traw, M] = size(Yraw);

% Generate lagged Y matrix
Ylag = mlag2(Yraw,p); % Y is [T x M]. ylag is [T x (Mp)]
X = [ones(Traw-p,1) Ylag(p+1:Traw,:)];
[Traw3, K] = size(X);
Y = Yraw(p+1:Traw,:); % This is the final Y matrix used for the VAR
T = Traw - p;

%--------------------------------PRIORS------------------------------------
% First get Ordinary Least Squares (OLS) estimators
A_OLS = inv(X'*X)*(X'*Y); % This is the matrix of regression coefficients
a_OLS = A_OLS(:);         % This is the vector of coefficients, i.e. it holds
                          % that a_OLS = vec(A_OLS)
SSE = (Y - X*A_OLS)'*(Y - X*A_OLS);
SIGMA_OLS = SSE./(T-K);

%-----------------Prior hyperparameters for bvar model
% Define hyperparameters
if prior == 1 % Noninformtive
    % I guess there is nothing to specify in this case!
    % Posteriors depend on OLS quantities
elseif prior == 2 % Minnesota
    A_prior = 0*ones(K,M);   
    a_prior = A_prior(:);
    
    % Hyperparameters on the Minnesota variance of alpha
    a_bar_1 = 1;
    a_bar_2 = 0.5;
    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./(T-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,:) = 1+i:M:K;
    end
   for i = 1:M  % for each i-th equation
        for j = 1:K   % for each j-th RHS variable
                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./(ceil((j-1)/M)^2); % variance on own lags           
                else
                    for kj=1:M
                        if find(j==ind(kj,:))>0
                            ll = kj;                   
                        end
                    end
                    V_i(j,i) = (a_bar_2*sigma_sq(i,1))./((ceil((j-1)/M)^2)*sigma_sq(ll,1));           
                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 a  
    
    % SIGMA is equal to the OLS quantity
    SIGMA = SIGMA_OLS;
    
elseif prior == 3 % Normal-Wishart (nat conj)
    % Hyperparameters on a ~ N(a_prior, SIGMA x V_prior)
    A_prior = 0*ones(K,M);   
    a_prior = A_prior(:);
    V_prior = 10*eye(K);
    % Hyperparameters on inv(SIGMA) ~ W(v_prior,inv(S_prior))
    v_prior = M;
    S_prior = eye(M);
    inv_S_prior = inv(S_prior);
end
    
%============================ POSTERIORS ==================================
%==========================================================================
    
%--------- Posterior hyperparameters of ALPHA and SIGMA with Diffuse Prior
if prior == 1
    % Posterior of alpha|Data ~ Multi-T(kron(SSE,inv(X'X)),alpha_OLS,T-K)
    V_post = inv(X'*X);
    a_post = a_OLS;
    A_post = reshape(a_post,K,M);
    
    % posterior of SIGMA|Data ~ inv-Wishart(SSE,T-K)
    S_post = SSE;
    v_post = T-K;
    
    % Now get the mean and variance of the Multi-t marginal posterior of alpha
    alpha_mean = a_post;
    alpha_var = (1/(v_post - M - 1))*kron(S_post,V_post);
    
%--------- Posterior hyperparameters of ALPHA and SIGMA with Minnesota Prior
elseif prior == 2
    % ******Get all the required quantities for the posteriors       
    V_post = inv( inv(V_prior) + kron(inv(SIGMA),X'*X) );
    a_post = V_post * ( inv(V_prior)*a_prior + kron(inv(SIGMA),X'*X)*a_OLS );
    A_post = reshape(a_post,K,M);
     
    % In this case, the mean is a_post and the variance is V_post
    alpha_mean=a_post;
%--------- Posterior hyperparameters of ALPHA and SIGMA with Normal-Wishart Prior
elseif prior == 3
    % ******Get all the required quantities for the posteriors       
    % For alpha
    V_post = inv( inv(V_prior) + X'*X );
    A_post = V_post * ( inv(V_prior)*A_prior + X'*X*A_OLS );
    a_post = A_post(:);
    
    % For SIGMA
    S_post = SSE + S_prior + A_OLS'*X'*X*A_OLS + A_prior'*inv(V_prior)*A_prior - A_post'*( inv(V_prior) + X'*X )*A_post;
    v_post = T + v_prior;
    
    % Now get the mean and variance of the Multi-t marginal posterior of alpha
    alpha_mean = a_post;
    alpha_var = (1/(v_post - M - 1))*kron(S_post,V_post); 
end

%======================= PREDICTIVE INFERENCE =============================
%==========================================================================

U = (Y - X*A_post)'; %The first two rows of U will be the reduced form residuals
SIGMA2 = U*U'/(T-1);% The upper left hand block gives the variance-covariance matrix

U = [U; zeros(M*(p-1),T)];
SIGMA = zeros(M*p,M*p);
SIGMA(1:M,1:M) = SIGMA2;
A = [A_post(2:end,:)'; [eye(M*(p-1)) zeros(M*(p-1),M)]];
V = [A_post(1,:)' ; zeros(M*p,1)]; %This gives the intercept terms



