分位数回归是计量，金融常用模型，无论是CoVaR还是VaR背后都有分位数回归的身影。我采用MCMC算法估计分位数回归模型，具体算法参考Korobilis(2017)文章。关于代码作以下说明：
1. 由于代码过长（完整代码接近800行），省略自己写的矩阵类（Matrix）和随机数类（Randgen）以及数据的载入、存储，头文件等，只粘贴估计过程（接近200行）
2.平时主要搞金融计量，matlab,C++通常只作为语言工具，因此代码疏漏，bug难免，欢迎C++高手批评指正
//////////////////////////////////////////////////////分隔符///////////////////////分隔符//////////////////////////////////////////////////////////
double **QuantileReg::Reg(double **Ydata, double **Xdata, double Prob, int T, int Vartotal,int Maxiter)
{
        Matrix M;//矩阵算法Class
        Randgen Gen;//随机数Class
        ////prior////
        double **Vinv, a0 = 0.1, b0 = 0.1;
        Vinv = new double*[Vartotal];
        for (int i = 0; i < Vartotal; i++)
        {
                Vinv[i] = new double[Vartotal];
                for (int j = 0; j < Vartotal; j++)
                {
                        if (i == j)
                        {
                                Vinv[i][j] = 1.0/9.0;
                        }
                        else
                        {
                                Vinv[i][j] = 0.0;
                        }
                }
        }
        /////preset//////
        double **beta, **z, sig2 = 0.1, theta = 0.0, tau = 0.0,tau_q,theta_q,*Err,**Tmp,b1,**Diag,sig2old=sig2;
        double **beta_post, **z_post, sig2_post = 0.0;
        int a1;
        beta = new double*[Vartotal];
        beta_post = new double*[Vartotal];
        for (int i = 0; i < Vartotal; i++)
        {
                beta[i] = new double[1];
                beta_post[i] = new double[1];
        }
        for (int i = 0; i < Vartotal; i++)
        {
                for (int j = 0; j < 1; j++)
                {
                        beta[i][j] = 0.0;
                        beta_post[i][j] = 0.0;
                }
        }
        z = new double*[T];
        z_post = new double*[T];
        for (int i = 0; i < T; i++)
        {
                z[i] = new double[1];
                z_post[i] = new double[1];
        }
        for (int i = 0; i < T; i++)
        {
                for (int j = 0; j < 1; j++)
                {
                        z[i][j] = 0.2;
                        z_post[i][j] = 0.0;
                }
        }
        Err = new double[T];
        Diag = new double*[T];
        for (int t = 0; t < T; t++)
        {
                Diag[t] = new double[T];
        }
        /////////////major cycle/////////////
        int maxiter = Maxiter, burnin = floor(maxiter*0.5);
        double **Xtile, **Ytile,**V1,**Mu1,**Rtmp;
        V1 = new double*[Vartotal];
        for (int i = 0; i < Vartotal; i++)
        {
                V1[i] = new double[Vartotal];
        }
        Ytile = new double*[T];
        for (int t = 0; t < T; t++)
        {
                Ytile[t] = new double[1];
        }
        for (int iter = 0; iter < maxiter; iter++)
        {
                std::cout << "迭代次数" << iter+1<< endl;
                 tau_q = 2.0 / (Prob*(1.0 - Prob));
                 theta_q = (1.0 - 2.0*Prob) / (Prob*(1.0 - Prob));
                 /////step1 drawing sig2///////
                 a1 = a0 + floor(1.5*T);
                 Tmp = M.Matrixdot(Xdata, beta, T, Vartotal, Vartotal, 1);
                 for (int t = 0; t < T; t++)
                 {
                         Err[t] = Ydata[t][0] - Tmp[t][0] - theta_q*z[t][0];
                         Err[t] = pow(Err[t], 2.0);
                 }
                 b1 = b0;
                 for (int t = 0; t < T; t++)
                 {
                         b1 = b1 + Err[t] / (2.0*z[t][0] + tau_q);
                 }
                 Tmp = Gen.IGrnd(a1, b1, 1, 1);
                 sig2 = Tmp[0][0];
                 if (sig2 < 1e-4)
                 {
                         sig2 = sig2old;
                 }
                 sig2old = sig2;
                 //////////step2 drawing beta///////////
                 for (int i = 0; i < T; i++)
                 {
                         for (int j = 0; j < T; j++)
                         {
                                 if (i == j)
                                 {
                                         Diag[i][j] = 1.0 / (sqrt(sig2)*tau_q*z[i][0]);
                                 }
                                 else
                                 {
                                         Diag[i][j] = 0.0;
                                 }
                         }
                 }
                 Tmp = M.Matrixtrans(Xdata, T, Vartotal);//矩阵转置
                 Xtile = M.Matrixdot(Tmp, Diag, Vartotal, T, T, T);//矩阵乘法
                 Tmp = M.Matrixdot(Xtile, Xdata, Vartotal, T, T, Vartotal);
                 for (int i = 0; i < Vartotal; i++)
                 {
                         for (int j = 0; j < Vartotal; j++)
                         {
                                 V1[i][j] = Tmp[i][j] + Vinv[i][j];
                         }
                 }
                 V1 = M.Matrixinv(V1, Vartotal);
                 if (_isnan(V1[Vartotal - 1][Vartotal - 1]))
                 {
                         V1 = M.Matrixinv(Vinv, Vartotal);
                 }
                 for (int t = 0; t < T; t++)
                 {
                         Ytile[t][0] = Ydata[t][0] - theta_q*z[t][0];
                 }
                 Tmp = M.Matrixdot(Xtile, Ytile, Vartotal, T, T, 1);
                 Mu1 = M.Matrixdot(V1, Tmp, Vartotal, Vartotal, Vartotal, 1);
                 Tmp = M.Matrixchol(V1, Vartotal);
                 Rtmp = Gen.Normrnd(Vartotal, 1);
                 Tmp = M.Matrixdot(Tmp, Rtmp, Vartotal, Vartotal, Vartotal, 1);
            for (int i = 0; i < Vartotal; i++)
                 {
                         beta[i][0] = Mu1[i][0] + Tmp[i][0];
                         //std::cout << beta[i][0] << "\t";
                 }
                 ////////////step3 draw z/////////////
                 Tmp = M.Matrixdot(Xdata,beta,T,Vartotal,Vartotal,1);
                 for (int t = 0; t < T; t++)
                 {
                         double k1 = sqrt(pow(theta_q, 2.0) + 2.0*tau_q) / abs(Ydata[t][0] - Tmp[t][0]);
                         double k2 = (pow(theta_q, 2.0) + 2.0*tau_q) / tau_q;
                         Rtmp = Gen.InvGrnd(k1, k2, 1, 1);
                         if (Rtmp[0][0]>1e-3)
                         {
                                 z[t][0] = Rtmp[0][0];
                         }
                         else
                         {
                                 z[t][0] = 1e-3;
                         }
                 }
                 /////////step4 save%%%%%%%%%%%%%%
                 if (iter >= burnin)
                 {
                         for (int i = 0; i < Vartotal; i++)
                         {
                                 beta_post[i][0] = beta_post[i][0] + beta[i][0] / (maxiter - burnin);
                         }
                         for (int t = 0; t < T; t++)
                         {
                                 z_post[t][0] = z_post[t][0] + z[t][0] / (maxiter - burnin);
                         }
                         sig2_post = sig2_post + sig2 / (maxiter - burnin);
                 }
        }
        ////////////释放内存///////////////////
        for (int t = 0; t < T; t++)
        {
                delete[] Diag[t];
        }
        delete[] Diag;
        for (int t = 0; t < Vartotal; t++)
        {
                delete[] Vinv[t];
                delete[] beta[t];
        }
        delete[] Vinv;
        delete[] beta;
        delete[] Err;
        for (int t = 0; t < T; t++)
        {
                delete[] Ytile[t];
                delete[] z[t];
        }
        delete[] Ytile;
        delete[] z;
        std::cout << "估计参数向量：" << endl;
        for (int i = 0; i < Vartotal; i++)
        {
                std::cout << beta_post[i][0] << "\t";
        }
        std::cout << "残差标准差：" << endl;
        std::cout << sig2_post << endl;
        return z;//返回潜变量序列
}
、