% MAIN_POSTERIOR0_CHOL.M 

clc;
clear all;
close all;
randn('seed',896)

load data_chol; [t,n]=size(y); p=12; H=17; T=t-p; 
M=2000;     %Number of reduced-form posterior draws

% Reduced form NIW prior and posterior
Y=y(p:t,:);	
for i=1:p-1
 	Y=[Y y(p-i:t-i,:)];
end
X=[ones(T,1) Y(1:t-p,:)];
Bhat=inv(X'*X)*X'*y(p+1:t,:);
B=Bhat';
Sigmahat = (y(p+1:t,:)-X*Bhat)'*(y(p+1:t,:)-X*Bhat)/T;
Ydep=y(p+1:t,:);

% Prior parameters
% Informative Minnesota prior for VAR model with intercept:
nu0=n+2;

% Coefficient prior centered on zero
Bbar0=zeros(p*n+1,n); 

% Residual variance of AR(1) for each variable
psi = NaN(n,1); 
for i=1:n
     ar1=ols1(y(2:end,i),[ones(t-1,1),y(1:end-1,i)]);
     psi(i)=ar1.sig2hatols;
end;
S0 = diag(psi);

% Prior coefficient variance
Omega=zeros(n*p+1,1); Omega(1)=1E6; d = n+2; lambda=0.2; alpha=2;
for i=1:p
    Omega(1+(i-1)*n+1:1+i*n)=(d-n-1)*(lambda^2)*(1/(i^alpha))./psi;
end
N0=diag(1./Omega);

% Posterior parameters
nuT   = T+nu0;
NT    = N0+X'*X;
BbarT = inv(NT)*(N0*Bbar0+X'*X*Bhat);
ST    = (nu0/nuT)*S0+(T/nuT)*Sigmahat+(1/nuT)*(Bhat-Bbar0)'*N0*inv(NT)*X'*X*(Bhat-Bbar0);
EvecB = reshape(BbarT,n*(n*p+1),1); 

Sigmamat = [];
for i=1:M
    RANTR=chol(inv(ST))'*randn(n,nuT)/sqrt(nuT);
    Sigma=inv(RANTR*RANTR');
    Sigmamat = [Sigmamat; reshape(Sigma,1,n*n)];    
end

% Shells
z=zeros(M,n+p*n^2+n^2+n^2);

% Evaluations of the posterior pdf at posterior draws
for j=1:M

    	disp(j)
    	Sigma  = reshape(Sigmamat(j,1:n*n),n,n);
    	VvecB = kron(Sigma,inv(NT));       
    	vecB   = EvecB+chol(VvecB)'*randn(n*(n*p+1),1);
    	A      = (reshape(vecB,1+n*p,n))'; 
    	P      = chol(Sigma)';     
      
        z(j,:) = [vecB' reshape(Sigma,1,n*n) reshape(P,1,n*n)];                     
         
end % end of j-loop

save posterior0_chol M z p n H t Ydep X;
