function [gammas,as,b0_tildes,acc,gs,sigmas,Qss,fore_results]=...
    MAIARSV(Y,qp,fblks,prior,c_metro,l_c,logfilepath,store_bvec,fore_opt)
% INPUTS:   Y is nx(T+p)      p is scalar (lags)

%           prior the structure containing
%           prior.gamma_bar, prior.Vgamma
%           prior.a_bar, prior.Va
%           prior.mb0,prior.vb0,
%           prior.g_bar, prior.Vg
%           prior.logsigma_bar ,prior.Ps_bar, prior.Qs_bar, prior.tau0s

%           l_c is the desired length of the chain
%           fblks is the 1xr  vector containing how many
%           components for each block [  n_1, n_2, ..., n_r ]  .
%           c_metro is a scaling factor to be calibrated so to have 30%
%           acceptance in the metropolis step
%           store_bvec is a 1x5 vector of booleans, to decide storage of elements:
%                   a, B0,  g, sigma, Q0s



%           fore_opt is the structure containing forecasting options (if
%           not provided the forest would not be computed). fore_opt.H and
%           foreopt.m are required variables, and Y_act as well

% Copyright (c) 2013-2021, Francesco Corsello, All Rights Reserved.
if nargin < 9
    bool_fore=false;
else
    bool_fore=true;
end

model=getenv('model');

fileID = fopen(logfilepath,'a');

%% data preparation
[mqp,iqp]=max(qp);
[y,~,~,~,m,n,T,x]=prep_data_TV_VAR(Y,mqp);  %y is nxT

%make_B0(b0_tilde,info_argument,model)

r=length(fblks);

q=qp(1); p=qp(2);
nrp=n*r*p;
rp=r*p;
nq=n*q;

x_bul=reshape(x(2:end,:),[n,mqp,T]);
x_bul_p=x_bul(:,1:p,:);
x_bul_q=x_bul(:,1:q,:);
inddiagXq=logical(repmat(eye(n),[1,q,T]));
calX_q=zeros(n,nq,T);
calX_q(inddiagXq)=x_bul_q(:);


%% MCMC settings
l_burn=round((3/7)*l_c);
l_tot=l_c+l_burn;
c_bar=1e-3;


%% pre allocation
%gamma, a, B0,g, sigma, Q0s
if store_bvec(1); gammas=zeros(nq,l_c);
else   gammas=[];   end
if store_bvec(2); as=zeros(nrp,l_c);
else   as=[];   end
if store_bvec(3); b0_tildes=zeros(size(prior.mb0,1),l_c);
else   b0_tildes=[];   end
if store_bvec(4); gs=zeros(m,l_c);
else   gs=[];   end
if store_bvec(5); sigmas=zeros(n,T,l_c);
else   sigmas=[];   end
if store_bvec(6); Qss=zeros(n,n,l_c);
else   Qss=[];   end

acc=zeros(size(prior.mb0,1),1);

%preallocate forecast (if requested)

if bool_fore
    fore_y_hat=zeros(n,fore_opt.H,l_c);
end

%% Initial conditions for the sampler

b_good_draw=0;

while not(b_good_draw)
    [~,a_h, ~,~, sigma_h , Qs]=draw_from_prior_MAIAR_SV(prior,T);
    b_good_draw=all(vec(sigma_h<50));
end
A_h=transpose(reshape(a_h,[rp,n]));
b0_tilde_h=prior.mb0;
g_h=0.01*randn(m,1);  %stable solution
G_h=make_tril_A_pages(g_h);
B0h=make_B0_restricted(b0_tilde_h,fblks);
Zh=reshape(mtimesx(B0h, x_bul_p),[rp,T]);    % rp xT

%%%%%%%%%%%%%%
%% some elements equal over the loop (efficiency gains)
precgamma=eye(nq)/prior.Vgamma;
wprecgammabar=precgamma*prior.gamma_bar;
iVA=eye(nrp)/prior.Va;
iVA_ma=iVA*prior.a_bar;
%%

disp('Gibbs sampler starts!')
%% Gibbs sampler
for ii=1:l_tot
    if mod(ii,1000)==0
        fprintf(fileID,['\n' datestr(now) ' Iteration number: ' num2str(ii)]);
    end
    %% ar coeffs, gammas
    y_circ=y-A_h*Zh;
    sigmainvrows_h=repmat(reshape(sigma_h.^(-1),[n,1,T]),[1,n,1]);
    SigmainvG_h=sigmainvrows_h.*repmat(G_h,[1,1,T]);  % n,n,T
    
    yctilSt=vec(mtimesx(SigmainvG_h,reshape(y_circ,[n,1,T]))); %nT x 1
    calXtilde_q=mtimesx(SigmainvG_h,calX_q);
    % calX_q   is n x nq x T
    calXSt=reshape(permute(calXtilde_q,[2,1,3]),[nq,n*T])'; %nT x nq
    XprX_=calXSt'*calXSt;
    precPost_=XprX_+precgamma;
    varPost_=precPost_\eye(nq);
    Xpry_=calXSt'*yctilSt;
    gamma_post=varPost_*(Xpry_+wprecgammabar);
    gamma_h=chol(varPost_,'lower')*randn(nq,1)+gamma_post;  %nq x 1
    y_noidio=y-reshape(mtimesx(calX_q,gamma_h),[n,T]);
    
    %%  loadings,   a
    G_inv_h=eye(n)/G_h;
    PAIh=CTA(y_noidio',Zh',n,rp,[],G_h,sigma_h',iVA,iVA_ma,A_h');
    a_h=vec(PAIh);
    %% b0_tilde
    A_h=transpose(reshape(a_h,[rp,n]));
    [b0_tilde_h,acc_h]=...
        RW_metro_step_b0('logLikKernel_MAI_VAR_SV',b0_tilde_h,fblks,y_noidio,A_h,x_bul_p,g_h,sigma_h,c_metro,prior.mb0,prior.vb0);
    B0h=make_B0_restricted(b0_tilde_h,fblks);
    Zh=reshape(mtimesx(B0h, x_bul_p),[rp,T]);
    
    %% off diagonal elements, G
    y_hat=y_noidio-A_h*Zh;
    W_h=create_offdiag_cov_regression(y_hat);
    
    Sigmasq_h=zeros(n,n,T);
    Sigmasq_h(logical(repmat(eye(n),[1,1,T])))=sigma_h(:).^2;
    
    g_h=sample_post_commcoeff_varspec_reg(y_hat,W_h, Sigmasq_h ,prior.g_bar, prior.Vg);
    G_h=make_tril_A_pages(g_h);
    
    %% sigmas
    y_tilde=G_h*y_hat;
    y_tilde_star=log( (y_tilde.^2) +c_bar);
    logsigma_h=log(sigma_h);
    
    [~,phi_h,Ups_h]=log_chisquared1_approx_states(y_tilde_star-2*logsigma_h);
    
    
    [logsigma_upd,logsigma_pred,Ps_upd,Ps_pred]=...
        FFBS_kfilter(y_tilde_star,phi_h,repmat(2*eye(n),[1,1,T]),Ups_h,...
        eye(n),Qs,prior.logsigma_bar,prior.Ps_bar);
    [logsigma_h,logsigma0]=...
        FFBS_sampling(logsigma_upd,logsigma_pred,Ps_upd,Ps_pred);
    
    sigma_h=exp(logsigma_h);
    
    %% sample innovation covariance matrix of stoch vols
    %sample from the IW
    Qs=draw_Q_from_SIW_posterior([logsigma0 logsigma_h], prior.Qs_bar, prior.tau0s);
    
    
    %% save steps after burn in
    
    if ii>l_burn
        ij=ii-l_burn;
        
        if store_bvec(1); gammas(:,ij)=gamma_h; end
        if store_bvec(2); as(:,ij)=a_h; end
        if store_bvec(3); b0_tildes(:,ij)=b0_tilde_h;  end
        if store_bvec(4); gs(:,ij)=g_h; end
        if store_bvec(5); sigmas(:,:,ij)=sigma_h; end
        if store_bvec(6); Qss(:,:,ij)=Qs; end
        
        acc=acc+acc_h./l_c;
        
        
        
        if bool_fore
            fore_y_hat(:,:,ij)=make_hstep_forecast_MAIAR(gamma_h,B0h,A_h,G_inv_h,sigma_h(:,end),y,qp,fore_opt.H);
        end
    end
    
    
    
end

if bool_fore
    fore_results.y_hat=fore_y_hat;
    fore_results.H=fore_opt.H;
else
    fore_results=struct([]);
end


fclose(fileID);
