function sample = multivariate_adjustment_gibbs(Y,X,sigma_alpha,ne,de,nw,dw,nv,dv,index_factors,nitermax);
% PURPOSE: Gibbs sampling multivariate adjustment model
% -----------------------------------------------------
% USAGE: 
% sample = multivariate_adjustment_gibbs(Y,X,sigma_alpha,ne,de,nw,dw,nv,dv,index_factors,nitermax)
% where
% sample = a structure containing the sample of the posterior distribution
% Y  (TxpxN) data dependent variables array with 
% T = size of the analysed period, 
% p = number of dependent variables
% N = number of firms
% X  (qxN) independent variables matrix with
% q = number of independent variables
% alphas(j) ~ N(0,sigma_alpha) j=1,...,N
% sigma_e(j) ~ IGamma(ne/2,de/2) j=1,...,N
% sigma_w(k,i) ~ IGamma(nw/2,dw/2) k=1,...,K (K = number of factors); i=1,...,N
% sigma_v(k,l) ~ IGamma(nv/2,dv/2) k,l=1,...,K;
% index_factors: binary matrix (pxK) that indicates, for each factor, the
% dependent variables related with it
% e.g. in the example analysed in the paper
% index_factors = 
%     1     0     0     0
%     1     0     0     0
%     1     0     0     0
%     0     1     0     0
%     0     1     0     0
%     0     0     1     0
%     0     0     1     0
%     0     0     0     1
% nitermax = maximum number of iterations
% -------------------------------------------------------------------------
% NOTE
% If there are not any independent variables use
% sample = multivariate_adjustment_gibbs(Y,X,sigma_alpha,ne,de,nw,dw,nv,dv,index_factors,nitermax);
% -------------------------------------------------------------------------
% RETURNS: a structure containing the sample of the posterior distribution
% sample.alpha = alphas;
% sample.beta = betas;
% sample.rho = rho;
% sample.phi = phi;
% sample.Ft = Ft;
% sample.gamma = gammas;
% sample.tau_e = taus_e;
% sample.taus_w = taus_w;
% sample.taus_v = taus_v;
%
% Initialising the output variables
%
N = size(Y,3);
p = size(Y,2);
T = size(Y,1);
K = size(index_factors,2);
q = size(X,1);
Z = zeros(q+1,N);
Z(1,:)=1;
if(q>0)
    Z(2:end,:)=X;
end
alphas = zeros(p,nitermax);
betas = ones(p,nitermax);
rho = zeros(K,K,nitermax);
phi = zeros(K,K,q,nitermax);
Ft = zeros(T+1,K,N);
gammas = zeros(K,K,N);
taus_e = ones(p,1);
taus_w = ones(K,N);
taus_v = ones(K,K);
index = 1:p;
factors = 1:K;
kfactors = zeros(p,1);
for(j=1:p)
    kfactors(j) = sum(index_factors(j,:).*factors);
end
variables_factors = sum(index_factors);
%
% Implementation of Gibbs sampling
%
for(it = 2:nitermax)
    it
    %
    % Update Ft
    % We use the forward filtering and backward sampling algorithm of Fruchwirth-Schnatter
    %
    %
    A1 = zeros(p,K);
    Sigma_e = diag(1./taus_e);
    for(k=1:K)
        index1 = index(index_factors(:,k)==1);
        A1(index1,k) =betas(index1,it-1);
    end
    for(i=1:N)
        % Forward filtering step
        A2 = reshape(gammas(:,:,i),K,K);
        mt = zeros(K,T+1);
        Ct = zeros(K,K,T+1);
        W = diag(1./taus_w(:,i));
        Rt = Ct;
        Ct(:,:,1) = eye(K);
        for(t=2:(T+1))
            at = A2*mt(:,t-1);
            Rt(:,:,t) = A2*reshape(Ct(:,:,t-1),K,K)*A2'+W;
            ft = A1*at;
            Qt = inv(A1*reshape(Rt(:,:,t),K,K)*A1'+Sigma_e);
            mt(:,t) = at+reshape(Rt(:,:,t),K,K)*A1'*Qt*(reshape(Y(t-1,:,i),p,1)-alphas(:,it-1)-ft);
            Ct(:,:,t) = reshape(Rt(:,:,t),K,K)-reshape(Rt(:,:,t),K,K)*A1'*Qt*A1*reshape(Rt(:,:,t),K,K);
        end
        % Backward sampling step
        aux = mvnrnd(mt(:,T+1),reshape(Ct(:,:,T+1),K,K))';
        Ft(T+1,:,i) = aux;
        for(t=T:(-1):1)
            var = reshape(Ct(:,:,t),K,K)*A2'*inv(reshape(Rt(:,:,t+1),K,K));
            med = (eye(K)-var*A2)*mt(:,t)+var*reshape(Ft(t+1,:,i),K,1);
            var = (eye(K)-var*A2)*reshape(Ct(:,:,t),K,K);
            aux = mvnrnd(med,var);
            Ft(t,:,i) = aux;
        end
    end
    %
    % Updating alpha, beta and tau_e
    %
    for(j=1:p)
        k = kfactors(j);
        if(variables_factors(k)>1)
           med = zeros(2,1);
            var = diag([1/sigma_alpha(j);1]);
            for(i=1:N)
                for(t=1:T)
                    med = med + [1;Ft(t+1,k,i)]*Y(t,j,i);
                    var = var + taus_e(j)*[1;Ft(t+1,k,i)]*[1;Ft(t+1,k,i)]';
                end
            end
            var = inv(var);
            med = var*(taus_e(j)*med);
            aux = mvnrnd(med,var);
            alphas(j,it) = aux(1);
            betas(j,it) = aux(2);
        else
            med = 0;
            var = 1/sigma_alpha(j) + N*T;
            for(i=1:N)
                for(t=1:T)
                    med = med + Y(t,j,i)-Ft(t+1,k,i);
                end
            end
            var = 1/var;
            med = var*(taus_e(j)*med);
            aux = normrnd(med,sqrt(var));
            alphas(j,it) = aux;
        end
        nt = ne+N*T;
        dt = de;
        for(i=1:N)
            for(t=1:T)
                dt = dt + (Y(t,j,i)-alphas(j,it)-betas(j,it)*Ft(t+1,k,i))^2;
            end
        end
        taus_e(j) = gamrnd(nt/2,2/dt);
    end
    % Updating gamma and tau_w 
    %
    Sigma_v = reshape(taus_v,K*K,1);
    Sigma_v = diag(Sigma_v);
    nt = nw + T-1;
    nt = nt*ones(K,1);
    for(i=1:N)
        Sigma_w = diag(taus_w(:,i));
        var = Sigma_v;
        med = reshape(rho(:,:,it-1),K*K,1);
        for(h=1:q)
            med = med + X(h,i)*reshape(phi(:,:,h,it-1),K*K,1);
        end
        med = Sigma_v*med;
        for(t=2:T)
            aux = reshape(Ft(t,:,i),K,1);
            var = var + kron(aux*aux',Sigma_w);
            aux1 = reshape(Ft(t+1,:,i),K,1);
            med = med + kron(aux,eye(K))*Sigma_w*aux1;
        end
        var = inv(var);
        med = var*med;
        aux = mvnrnd(med,var);
        gammas(:,:,i) = reshape(aux,K,K);
        dt = dw*ones(K,1);
        for(t=2:T)
            dt = dt + (reshape(Ft(t+1,:,i),K,1)-reshape(gammas(:,:,i),K,K)*reshape(Ft(t,:,i),K,1)).^2;
        end
        taus_w(:,i) = gamrnd(nt/2,2./dt); 
    end
    %
    % Updating rho and phi
    %
    var = Z*Z';
    med = zeros((q+1)*K*K,1);
    for(i=1:N)
        g = reshape(gammas(:,:,i),K*K,1);
        med = med + kron(Z(:,i),eye(K*K))*Sigma_v*g;
    end
    var = eye((q+1)*K*K)+kron(var,Sigma_v);
    var = inv(var);
    med = var*med;
    aux = mvnrnd(med,var);
    aux = reshape(aux,K*K,q+1);
    rho(:,:,it) = reshape(aux(:,1),K,K);
    for(j=2:(q+1))
        phi(:,:,j-1,it) =reshape(aux(:,j),K,K);
    end
    % Updating tau_v
    %
    nt = nv+N;
    for(j=1:K)
        for(k=1:K)
            dt=dv;
            for(i=1:N)
                if(q>0)
                    dt = dt + (gammas(j,k,i)-rho(j,k,it)-reshape(phi(j,k,:,it),1,q)*X(:,i)).^2;
                else
                    dt = dt + (gammas(j,k,i)-rho(j,k,it)).^2;
                end
            end
            taus_v(j,k) = gamrnd(nt/2,2/dt);
        end
    end
end
%
% Results
%
sample.alpha = alphas;
sample.beta = betas;
sample.rho = rho;
sample.phi = phi;
sample.Ft = Ft;
sample.gamma = gammas;
sample.tau_e = taus_e;
sample.taus_w = taus_w;
sample.taus_v = taus_v;

    
    
                