function [beta_keep, Sigma_keep]=Multi_fit(y, x, nrep, warmup, T, k, num_eq)

yt=y';         % Yt:  num_eq x T           
for i=1:T       % Xt: num_eq x k x T, three dimension matrix
    xt(:,:,i)=kron(eye(num_eq), x(i,:)); 
end

% First, set some priors for coefficients and the covariance matrix b~N(ub,vb), H^-1 ~W(omeg, v)
ubeta=zeros(num_eq*k,1); vbeta=5*eye(num_eq*k);
Sigma_uncond=cov(y);
v0=20; omeg0=inv(v0.*Sigma_uncond);     % Sigma0_draw ~ Inv(Wishrnd(omeg0, v0))
beta_draw=multi_normrnd(ubeta, vbeta,1)';
%Sigma_draw=inv(wishrnd(omeg0,v0));
Sigma_draw=Sigma_uncond;
%% Gibbs-sampling
for iter=1:nrep+warmup
    % Step 1: sample the beta from normal distribution
    Db=zeros(num_eq*k); db=zeros(num_eq*k,1);
    for i=1:T
        Db=Db+ xt(:,:,i)'*inv(Sigma_draw)*xt(:,:,i);
        db=db+xt(:,:,i)'*inv(Sigma_draw)*yt(:,i);
    end    
    Db=(Db+inv(vbeta))^-1;
    db=db+inv(vbeta)*ubeta;
    beta_draw=multi_normrnd(Db*db, Db,1)';
 
    % sample the inverse Sigma matrix~ Wishart distribution~W(post_omeg, post_v)
    post_omeg=zeros(num_eq);
    for i=1:T
        post_omeg=post_omeg+(yt(:,i)- xt(:,:,i)*beta_draw)*(yt(:,i)- xt(:,:,i)*beta_draw)';
    end    
    post_omeg=inv(post_omeg+inv(omeg0));
    post_v=T+v0;
    Sigma_draw= inv(wishrnd(post_omeg,post_v));

    if iter>warmup
        beta_keep(iter-warmup,:)=beta_draw';
        Sigma_keep(:,:, iter-warmup)=Sigma_draw;
    end
    fprintf(' Iteration to %3g \n', iter)
end















