function Samples = mcmc_factor(Factors_Return,Opts)
%**************************************************************************
% MATLAB/Octave用プログラム（MATLABではStatistics Toolboxが必要）
% ギブズ・サンプラーによるファクターの分布のパラメータの
% モンテカルロ標本の生成
% 慶應義塾大学経済学部　中妻照雄
%**************************************************************************
% T --- 期関数
% N --- 銘柄数
% K --- ファクター数
% 入力
% Factors_Return --- ファクター (T x K)
% Opts --- ハイパーパラメータなどの設定を含む構造体
% (1) Factor_Mean_Mean ---  ファクターの期待値の平均 (K x 1)
% (2) Factor_Mean_Var --- ファクターの期待値の分散共分散行列 (K x K)
% (3) Factor_Var_Shape --- ファクターの分散共分散行列の形状パラメータ
% (4) Factor_Var_Scale --- ファクターの分散共分散行列の規模行列 (K x K)
% (5) Burnin ---  検査稼働の回数
% (6) Sample --- ギブズ・サンプラーの繰り返し回数
% (7) Batch --- サンプルの抽出間隔
% 出力
% Samples --- モンテカルロ標本の構造体
% (1) Factor_Mean --- ファクターの期待値のモンテカルロ標本
% (2) Factor_Var --- ファクターの標準偏差と相関係数のモンテカルロ標本
%**************************************************************************
[T,K] = size(Factors_Return);
vech_idx_K = reshape(tril(reshape((1:K*K),K,K)),1,K*K);
vech_idx_K = vech_idx_K(vech_idx_K>0);
IOmega_F = Opts.Factor_Mean_Var\eye(K);
IOmegaMu_F = IOmega_F*Opts.Factor_Mean_Mean;
Nu_F = Opts.Factor_Var_Shape+T;
DF_F = .5*linspace(Nu_F,Nu_F-K+1,K);
Scale_F = Opts.Factor_Var_Scale;
F = Factors_Return';
sum_F = sum(F,2);
bar_F = mean(F,2);
% ギブズ・サンプラー
Sample_Size = Opts.Sample;
Burnin_Size = Opts.Burnin;
Batch_Size = Opts.Batch;
Iterations = Sample_Size+Burnin_Size;
Factor_Mean = zeros(floor(Sample_Size/Batch_Size),K);
Factor_Var = zeros(floor(Sample_Size/Batch_Size),K*(K+1)/2);
jt = 1;
for it = 1:Iterations
    % ファクターの分散共分散行列$\Sigma_F$の生成
    U = F-repmat(bar_F,1,T);
    S_F = Scale_F+U*U';
    L_F = chol(S_F\eye(K))';
    A_F = tril(randn(K,K),-1)+diag(sqrt(2*randg(DF_F)));
    ISigma_F = L_F*(A_F*A_F')*L_F';
    % MATLABのStatistics Toolboxでは関数wishrnd()を利用できる
    % その場合は49-51行を削除して下の行に置き換える
    %ISigma_F = wishrnd(S_F\eye(K),Nu_F);
    % ファクターの平均$\lambda$の生成
    bar_F_Var = (T*ISigma_F+IOmega_F)\eye(K);
    bar_F_Mean = bar_F_Var*(ISigma_F*sum_F+IOmegaMu_F);
    bar_F = bar_F_Mean+chol(bar_F_Var)'*randn(K,1);
    % モンテカルロ標本の保存
    if it > Burnin_Size && mod(it,Batch_Size)==0
      Factor_Mean(jt,:) = bar_F;
      Sigma_F = ISigma_F\eye(K);
      SD_F = sqrt(diag(Sigma_F));
      Corr_F = diag(1./SD_F)*Sigma_F*diag(1./SD_F);
      Cov_F = tril(Corr_F,-1)+diag(SD_F);
      Factor_Var(jt,:) = Cov_F(vech_idx_K);
      jt = jt+1;
    end
end
Samples.Factor_Mean = Factor_Mean;
Samples.Factor_Var = Factor_Var;
