function Samples = mcmc_mfm(Stocks_Return,Factors_Return,Opts)
%**************************************************************************
% MATLAB/Octave用プログラム（MATLABではStatistics Toolboxが必要）
% ギブズ・サンプラーによるマルチファクター・モデルのパラメータの
% モンテカルロ標本の生成
% 慶應義塾大学経済学部　中妻照雄
%**************************************************************************
% T --- 期関数
% N --- 銘柄数
% K --- ファクター数
% 入力
% Stocks_Return  --- 個別銘柄の収益率 (T x N)
% Factors_Return --- ファクター (T x K)
% Opts --- ハイパーパラメータなどの設定を含む構造体
% (1) Coef_Mean_Mean ---  回帰係数の期待値の平均 (K+1 x 1)
% (2) Coef_Mean_Var --- 回帰係数の期待値の分散共分散行列 (K+1 x K+1)
% (3) Coef_Var_Shape --- 回帰係数の分散共分散行列の形状パラメータ
% (4) Coef_Var_Scale --- 回帰係数の分散共分散行列の規模行列 (K+1 x K+1)
% (5) Error_Var_Shape --- 誤差項の分散共分散行列の形状パラメータ
% (6) Error_Var_Scale --- 誤差項の分散共分散行列の規模行列 (N x N)
% (7) Leverage --- 自己資本比率 (1 x N)
% (8) Burnin ---  検査稼働期間の回数
% (9) Sample --- ギブズ・サンプラーの繰り返し回数
% (10) Batch --- サンプルの抽出間隔
% 出力
% Samples --- モンテカルロ標本の構造体
% (1) MFM_Coef --- 回帰係数のモンテカルロ標本
%                  最初のK列は業種平均回帰係数
%                  それ以降は第1企業，第2企業，．．．の回帰係数
% (2) MFM_Var --- 分散共分散行列のモンテカルロ標本
%                 最初の(K+1)(K+2)/2列は回帰係数の標準偏差と相関係数
%                 残りのN(N+1)/2列は誤差項の標準偏差と相関係数
%**************************************************************************
[T,N] = size(Stocks_Return); 
K = size(Factors_Return,2)+1; % Kは回帰係数の数 (ファクター数 + 1)
KN = K*N;
vech_idx_K = reshape(tril(reshape((1:K*K),K,K)),1,K*K);
vech_idx_K = vech_idx_K(vech_idx_K>0);
vech_idx_N = reshape(tril(reshape((1:N*N),N,N)),1,N*N);
vech_idx_N = vech_idx_N(vech_idx_N>0);
IOmega_B = Opts.Coef_Mean_Var\eye(K);
IOmegaMu_B = IOmega_B*Opts.Coef_Mean_Mean;
Nu_E = Opts.Error_Var_Shape+T;
DF_E = .5*linspace(Nu_E,Nu_E-N+1,N);
Nu_B = Opts.Coef_Var_Shape+N;
DF_B = .5*linspace(Nu_B,Nu_B-K+1,K);
Scale_E = Opts.Error_Var_Scale;
Scale_B = Opts.Coef_Var_Scale;
Lev = [ones(1,N); repmat(Opts.Leverage,K-1,1)];
Y = Stocks_Return;
X = [ones(T,1) Factors_Return];
XX = X'*X;
XY = X'*Y;
B = X\Y;
bar_B = mean(B,2);
U = B-repmat(bar_B,1,N);
ISigma_B = (Nu_B-K-1)*((U*U'+Scale_B)\eye(K));
% ギブズ・サンプラー
Sample_Size = Opts.Sample;
Burnin_Size = Opts.Burnin;
Batch_Size = Opts.Batch;
Iterations = Sample_Size+Burnin_Size;
MFM_Coef = zeros(floor(Sample_Size/Batch_Size),K+KN);
MFM_Var = zeros(floor(Sample_Size/Batch_Size),K*(K+1)/2+N*(N+1)/2);
jt = 1;
for it = 1:Iterations
    % 誤差項の分散共分散行列$\Sigma_E$の生成
    E = Y-X*B;
    S_E = Scale_E+E'*E;
    L_E = chol(S_E\eye(N))';
    A_E = tril(randn(N,N),-1)+diag(sqrt(2*randg(DF_E)));
    ISigma_E = L_E*(A_E*A_E')*L_E';
    % MATLABのStatistics Toolboxでは関数wishrnd()を利用できる
    % その場合は70-72行を削除して下の行に置き換える
    %ISigma_E = wishrnd(S_E\eye(N),Nu_E);
    % 回帰係数$B$の生成
    B_Var = (kron(ISigma_E,XX)+diag(reshape(Lev,KN,1)) ...
	     *kron(eye(N),ISigma_B)*diag(reshape(Lev,KN,1)))\eye(KN);
    B_Mean = B_Var*reshape((XY*ISigma_E+Lev.*repmat(ISigma_B*bar_B,1,N)),KN,1);
    b_v = B_Mean+chol(B_Var)'*randn(KN,1);
    B = reshape(b_v,K,N);
    LevB = Lev.*B;
    % 回帰係数$B$の分散共分散行列$\Sigma_B$の生成
    U = LevB-repmat(bar_B,1,N);
    S_B = Scale_B+U*U';
    L_B = chol(S_B\eye(K))';
    A_B = tril(randn(K,K),-1)+diag(sqrt(2*randg(DF_B)));
    ISigma_B = L_B*(A_B*A_B')*L_B';
    % MATLABのStatistics Toolboxでは関数wishrnd()を利用できる
    % その場合は86-88行を削除して下の行に置き換える
    %ISigma_B = wishrnd(S_B\eye(K),Nu_B);
    % 回帰係数$B$の平均$\bar B$の生成
    bar_B_Var = (N*ISigma_B+IOmega_B)\eye(K);
    bar_B_Mean = bar_B_Var*(ISigma_B*sum(LevB,2)+IOmegaMu_B);
    bar_B = bar_B_Mean+chol(bar_B_Var)'*randn(K,1);
    % モンテカルロ標本の保存
    if it > Burnin_Size && mod(it,Batch_Size)==0
      MFM_Coef(jt,:) = [bar_B' reshape(B,1,KN)];
      Sigma_B = ISigma_B\eye(K);
      SD_B = sqrt(diag(Sigma_B));
      Corr_B = diag(1./SD_B)*Sigma_B*diag(1./SD_B);
      Cov_B = tril(Corr_B,-1)+diag(SD_B);
      vech_B = Cov_B(vech_idx_K);
      Sigma_E = ISigma_E\eye(N);
      SD_E = sqrt(diag(Sigma_E));
      Corr_E = diag(1./SD_E)*Sigma_E*diag(1./SD_E);
      Cov_E = tril(Corr_E,-1)+diag(SD_E);
      vech_E = Cov_E(vech_idx_N);
      MFM_Var(jt,:) = [vech_B vech_E];
      jt = jt+1;
    end
end
Samples.MFM_Coef = MFM_Coef;
Samples.MFM_Var = MFM_Var;
