精华内容
下载资源
问答
  • matlab VAR模型应用实例,附源代码和PPT
  • TVP-VAR模型MATLAB代码

    2019-05-18 11:05:20
    TVP-VAR模型MATLAB代码,修改变量与数据可以直接运行,很方便!
  • matlab中的MVAR模型代码S-MVAR Matlab 工具箱允许分析计算 VAR 模型的参数,探索稀疏回归和状态空间 (SS) 模型的组合方法。 具体而言,研究的方法是:普通最小二乘分析、LASSO 回归、弹性网络回归、融合 LASSO 回归...
  • MATLABvar模型代码VAE-GAN-CelebA 与本文相关的Python代码:VanRullen&Reddy(2019) 该文件夹包含: 一个(在CelebA数据集上约15个时期= 50,000个批次的64张图像) VAE-GAN人脸分解/重建模型的一组.py函数,尤其...
  • 金融风险VaR模型研究\蒙特卡罗算法与matlab_精品教程_.pdf
  • 金融风险度量的 VaR 模型 作 者 郝芸芸 系 别 经济管理学院 专 业 金融学 学 号 31403020 摘要VaR 是使投资风险数量化的工具旨在估计给定金融资产 组合在正常的资产价格波动下未来可能的或潜在的损失目前常用 的 VaR...
  • MATLABvar模型代码二阶信息有助于大规模的视觉识别吗? 由和创建 内容 介绍 该存储库包含在ImageNet 2012数据集上针对以下论文训练的源代码和模型: @article{Li2017, author = {Peihua Li,Jiangtao Xie,Qilong ...
  • 强调 : 从 FRED 加载数据并转换数据以获得平稳性将转换后的数据划分为预采样、估计和预测区间制作多个模型模型拟合到数据使用各种回测技术确定模型根据最佳模型进行预测 产品重点: 的MATLAB DataFeed 工具箱...
  • TVP-VAR模型的贝叶斯估计 此仓库包含有关如何使用TVP-VAR模型进行贝叶斯分析的信息。 在深入研究代码之前,您应该看一下Bayes_TVPVAR_Presentation文件。 这将使您对TVP-VAR与常规VAR模型有何不同以及我们如何在TVP...
  • cvar代码matlab 多尺度,多目标星模型检测算法 多尺度多目标星形模型Matlab代码,用于单发目标检测和识别。 每个对象使用一个或几个示例进行训练,一次即可支持多达数千个对象。 作者:Leonid Karlinsky(),Joseph...
  • 一、简介 VAR:向量自回归模型,结果仅具有统计上的 意义 SVAR:结构向量自回归模型 TVP-VAR:Time Varying Parameter-Stochastic Volatility-Vector Auto Regression。...在VAR模型中,所有的参数遵循一阶游走

    一、简介

    VAR:向量自回归模型,结果仅具有统计上的 意义
    SVAR:结构向量自回归模型
    TVP-VAR:Time Varying Parameter-Stochastic Volatility-Vector Auto Regression。时变参数随机波动率向量自回归模型,与VAR 不同的是,模型没有同方差的假定,更符合实际。并且时变参数假定随机波动率,更能捕捉到经济变量在不同时代背景下所具有的关系和特征(时变影响)。将随机波动性纳入TVP估算中可以显着提高估算性能。在VAR模型中,所有的参数遵循一阶游走过程。
    随机波动的概念在TVP-VAR中很重要,随机波动是1976年由Black提出,随后在经济计量中有很大的发展。近几年,随机波动也经常被用在宏观经济的经验分析中。很多情况下,经济数据的产生过程中具有漂移系数和随机波动的冲击,如果是这种情况,那么使用具有时变系数但具有恒定波动性的模型会引起一个问题,即由于忽略了扰动中波动性的可能变化,估计的时变系数可能会出现偏差。为了避免这种问题,TVP-VAR模型假设了随机波动性。尽管似然函数难以处理,所以随机波动率使估计变得困难,但是可以在贝叶斯推断中使用马尔可夫链蒙特卡罗(MCMC)方法来估计模型。
    从 具有随机波动性的时变参数(TVP)回归模型的估计算法(TVP-VAR模型的单变量情况)说起。
    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述
    对于MCMC算法,π(θ)为先验密度,后验分布 π(θ,α,h|y)^4。
    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述
    四、进行TVP-VAR建模时,也需要数据平稳。可以用ADF单位根检验法检验数据的平稳性,不平稳的数据可以做差分,直至平稳,用差分后的数据进行建模。

    需要注意的是:TVP-VAR建模时,变量之间的顺序会影响到最后的实证结果。变量的顺序问题可能是实证结果不符合预期的一个原因。最好把关注的变量放在首位,这样在时变关系图中就能得到较好的展示。

    二、源代码

    %==================================4.1 Main================================
    clear;
    clc;
     
    % Load Korobilis (2008) quarterly data
    load ydata.dat; % data 
    load yearlab.dat; % data labels 
     
    %%
    %----------------------------------BASICS----------------------------------
     
     
    Y=ydata;
    t=size(Y,1);        % t - The total number of periods in the raw data (t=215)
    M=size(Y,2);        % M - The dimensionality of Y (i.e. the number of variables)(M=3)
    tau = 40;           % tau - the size of the training sample (the first forty quarters)
    p = 2;              % p - number of lags in the VAR model 
     
    %% Generate the Z_t matrix, i.e. the regressors in the model. 
     
    ylag = mlag2(Y,p);          % This function generates a 215x6 matrix with p lags of variable Y. 
    ylag = ylag(p+tau+1:t,:);   % Then remove our training sample, so now a 173x6 matrix. 
     
    K = M + p*(M^2);            % K is the number of elements in the state vector
     
    % Here we distribute the lagged y data into the Z matrix so it is
    % conformable with a beta_t matrix of coefficients. 
     
    Z = zeros((t-tau-p)*M,K);   
    for i = 1:t-tau-p
        ztemp = eye(M);
        for j = 1:p        
            xtemp = ylag(i,(j-1)*M+1:j*M); 
            xtemp = kron(eye(M),xtemp);
            ztemp = [ztemp xtemp];
        end
        Z((i-1)*M+1:i*M,:) = ztemp;
    end
     
    % Redefine our variables to exclude the training sample and the first two
    % lags that we take as given, taking total number of periods (t) from 215 
    % to 173. 
     
    y = Y(tau+p+1:t,:)';
    yearlab = yearlab(tau+p+1:t);
    t=size(y,2); % t now equals 173
     
    %% --------------------MODEL AND GIBBS PRELIMINARIES-----------------------
     
    nrep = 5000;    % Number of sample draws 
    nburn = 2000;   % Number of burn-in-draws
    it_print = 100; % Print in the screen every "it_print"-th iteration
     
     
    %% INITIAL STATE VECTOR PRIOR  
     
    % We use the first 40 observations (tau) to run a standard OLS of the
    % measurement equation, using the function ts_prior. The result is
    % estimates for priors for B_0 and Var(B_0). 
     
    [B_OLS,VB_OLS]= ts_prior(Y,tau,M,p);
     
    % Given the distributions we have, we now have to define our priors for B,
    % Q and Sigma. These are set in accordance with how they are set in 
    % Primiceri (2005). These are the hyperparameters of the beta, Q and Sigma
    % initial priors. 
     
    B_0_prmean = B_OLS;
    B_0_prvar = 4*VB_OLS; 
     
    Q_prmean = ((0.01).^2)*tau*VB_OLS;
    Q_prvar = tau;
     
    Sigma_prmean = eye(M);
    Sigma_prvar = M+1;
     
    % To start the Kalman filtering assign arbitrary values that are in support
    % of their priors, Q and Sigma. 
     
    consQ = 0.0001;
    Qdraw = consQ*eye(K);
    Sigmadraw = 0.1*eye(M);
     
    % Create some matrices for storage that will be filled in once we
    % start the Gibbs sampling.
     
    Btdraw = zeros(K,t); 
    Bt_postmean = zeros(K,t);
    Qmean = zeros(K,K);
    Sigmamean = zeros(M,M);
     
    %% -------------------------IRF-PRELIMINARIES------------------------------
    nhor = 21;     % The number of periods in the impulse response function. 
     
    % Matricies to be filled containing IRFs for 1975q1, 1981q3, 1996q1. The 
    % dimensions correspond to the iterations of the gibbs sample, each of the
    % variables, and each of the 21 periods of the IRF analysis. 
     
    imp75 = zeros(nrep,M,nhor); 
    imp81 = zeros(nrep,M,nhor);
    imp96 = zeros(nrep,M,nhor);
     
    % This corresponds to variable J introduced in equation (14) in the report
     
    bigj = zeros(M,M*p); 
    bigj(1:M,1:M) = eye(M);
     
     
    %% ================ START GIBBS SAMPLING ==================================
     
    tic; % This is just a timer
    disp('Number of iterations');
     
     
    for irep = 1:nrep + nburn    % 7000 gibbs iterations starts here
        % Print iterations - this just updates on the progress of the sampling
        if mod(irep,it_print) == 0
            disp(irep);toc;
        end
        
    %% Draw 1: B_t from p(B_t|y,Sigma)
        
        % We use the function 'carter_kohn_hom' to to run the FFBS algorithm. 
        % This results in a 21x173 matrix, corresponding to one Gibbs sample
        % draw of each of the coefficients in each time period. The inputs
        % Sigmadraw and Qdraw are updated for each Gibbs sample repetition.
        
        [Btdraw] = carter_kohn_hom(y,Z,Sigmadraw,Qdraw,K,M,t,B_0_prmean,B_0_prvar);
        
      
    %% Draw 2: Q from p(Q^{-1}|y,B_t) which is i-Wishart
     
        % We draw Q from an Inverse Wishart distribution. The parameters 
        % of the distribution are derived as equation (11) in the main report.
        % The mean is taken as the inverse of the accumulated sum of squared 
        % errors added to the prior mean, and the variance is simply t.  
        
        % Differencing Btdraw to create the sum of squared errors
        
        Btemp = Btdraw(:,2:t)' - Btdraw(:,1:t-1)'; 
        sse_2Q = zeros(K,K);
        for i = 1:t-1
            sse_2Q = sse_2Q + Btemp(i,:)'*Btemp(i,:);
        end
     
        Qinv = inv(sse_2Q + Q_prmean);      % compute mean to use for Wishart draw
        Qinvdraw = wish(Qinv,t+Q_prvar);    % draw inv q from the wishart distribution
        Qdraw = inv(Qinvdraw);              % find non-inverse q 
        
    %% Draw 3: Sigma from p(Sigma|y,B_t) which is i-Wishart
     
        % We draw Sigma from an Inverse Wishart distribution. The parameters 
        % of the distirbution are derived as equation (10) in the main report. 
        % The mean is taken as the inverse of the sum of squared residuals
        % added to the prior mean. The variance is simply t. 
        
        % Find residuals using data and the current draw of coefficients
        
        resids = zeros(M,t);
        for i = 1:t
            resids(:,i) = y(:,i) - Z((i-1)*M+1:i*M,:)*Btdraw(:,i);
        end
        
        % Create a matrix for the accumulated sum of squared residuals, to
        % be used as the mean parameter in the i-Wishart draw below. 
        sse_2S = zeros(M,M);
        for i = 1:t
            sse_2S = sse_2S + resids(:,i)*resids(:,i)';
        end
        
        Sigmainv = inv(sse_2S + Sigma_prmean);          % compute mean to use for the Wishart
        Sigmainvdraw = wish(Sigmainv,t+Sigma_prvar);    % draw from the Wishsart distribution
        Sigmadraw = inv(Sigmainvdraw);                  % turn into non-inverse Sigma
        
    %% IRF 
        % We only apply IRF analysis once we have exceeded the burn-in draws phase.
        if irep > nburn;         
                
                % Create matrix that is going to contain all beta draws over
                % which we will take the mean after the Gibbs sampler as our moment
                % estimate: 
                Bt_postmean = Bt_postmean + Btdraw;
                
                % biga is the A matrix of the VAR(1) version of our VAR(2) model, 
                % found in equation (12. biga changes in every period of the 
                % analysis, because the coefficients are time varying, so we 
                % apply the analysis below in every time period. 
                
                biga = zeros(M*p,M*p); 
                for j = 1:p-1
                    biga(j*M+1:M*(j+1),M*(j-1)+1:j*M) = eye(M); % fill the A matrix with identity matrix (3) in bottom left corner
                end
     
                % The following procedure is applied separately in each time period. 
                
                % This loop takes coefficients of the relevant time period from
                % Bt_draw (which contains all coefficients for all t) and uses
                % them to update the biga matrix, so that it can change for
                % every t. 
                
                for i = 1:t 
                    bbtemp = Btdraw(M+1:K,i);  % get the draw of B(t) at time i=1,...,T  (exclude intercept)
                    splace = 0;
                    for ii = 1:p
                        for iii = 1:M
                            biga(iii,(ii-1)*M+1:ii*M) = bbtemp(splace+1:splace+M,1)'; % Load non-intercept coefficient draws
                            splace = splace + M;
                        end
                    end
                    
                    % Next we want to create a shock matrix in which the third
                    % column is [0 0 1]', therefore implementing a unit shock
                    % in the interest rate. 
                    
                    shock = eye(3);
                    
                    % Now get impulse responses for 1 through nhor future
                    % periods. impresp is a 3x63 matrix which contains 9
                    % response values in total for each period, 3 for each 
                    % variable. These three responses correspond to the 3
                    % possible shocks that are contained in the schock
                    % matrix
                    
                    % bigai is updated through mulitiplication with the 
                    % coefficient matrix after each time period. 
                                    
                    % Create a results matrix to store impulse responses in all periods
                    impresp = zeros(M,M*nhor); 
                    
                    % Fill in the first period of the results matrix with the shock (as defined above) 
                    impresp(1:M,1:M) = shock;
                    
                    % Create a separate variable for the a matrix so that we
                    % can update it for each period of the IRF analysis. 
                    bigai = biga; 
                    
                    % This follows the impulse response function as in equation 15.
                    % Fill in each period of the results matrix according to
                    % the impulse response function formula. 
                    
                    for j = 1:nhor-1
                        impresp(:,j*M+1:(j+1)*M) = bigj*bigai*bigj'*shock;
                        bigai = bigai*biga; % update the coefficient matrix for next period
                    end
     
                    % The section below keeps only the responses that we are interested in:
                    % - those from the periods 1975q1, 1981q3, and 1996q1
                    % - those that correspond to the shock in the interest 
                    % rate (i.e. those caused by the third column of our shock
                    % matrix). 
                    
                    if yearlab(i,1) == 1975.00;   % store only IRF from 1975:Q1
                        impf_m = zeros(M,nhor);
                        jj=0;
                        for ij = 1:nhor
                            jj = jj + M;    % select only the third column for each time period of the IRF
                            impf_m(:,ij) = impresp(:,jj);
                        end           
                    % For each iteration of the Gibbs sampler, fill in the
                    % results along the first dimension 
                        imp75(irep-nburn,:,:) = impf_m; 
                    end
                    
                    if yearlab(i,1) == 1981.50;   % store only IRF from 1981:Q3
                        impf_m = zeros(M,nhor);
                        jj=0;
                        for ij = 1:nhor
                            jj = jj + M;    % select only the third column for each time period of the IRF
                            impf_m(:,ij) = impresp(:,jj);
                        end
                    % For each iteration of the Gibbs sample, fill in the
                    % results along the first dimension 
                        imp81(irep-nburn,:,:) = impf_m; 
                    end
                    
                    if yearlab(i,1) == 1996.00;   % store only IRF from 1996:Q1
                        impf_m = zeros(M,nhor);
                        jj=0;
                        for ij = 1:nhor
                            jj = jj + M;    % select only the third column for each time period of the IRF
                            impf_m(:,ij) = impresp(:,jj);
                        end
                    % For each iteration of the Gibbs sample, fill in the
                    % results along the first dimension 
                        imp96(irep-nburn,:,:) = impf_m; 
                    end
                end % End getting impulses for each time period 
            end % End the impulse response calculation section   
    end % End main Gibbs loop (for irep = 1:nrep+nburn)
     
    clc;
    toc; % Stop timer and print total time
    %% ================ END GIBBS SAMPLING ==================================
     
    % Even though it is not used in our IRF analysis since we are integrating
    % that into the Gibbs sampling loop, here is how to take the mean of the 
    % draw of the betas as moment estimate:
     
    Bt_postmean = Bt_postmean./nrep;
     
    %% Graphs and tables
     
    % This works out the percentage range of that each variable's coefficient
    % spans across time. I.e. to what extent was the TVP facility used by each
    % variable in the model? This is calculated by finding the range for each
    % variable as a percentage of the mean coefficient size for that variable.
    % The result is a 21x1 vector, and it is reshaped into a 3x7 matrix in order
    % to map onto the system of equations (2,3,and 4) set out in the report.
     
    Bt_range = ones(21,1)
    for i = 1:21
        Bt_range(i) = abs((max(Bt_postmean(i,:))-min(Bt_postmean(i,:)))/mean(Bt_postmean(i,:)))
    end 
    Bt_range = reshape(Bt_range,3,7)
     
    % Create a table of coefficient ranges for export to the report
     
    rowNames = {'Inflation','Unemployment','Interest Rate'};
    colNames = {'Intercept','Inf_1','Unemp_1', 'IR_1','Inf_2','Unemp_2', 'IR_2'};
    pc_change_table = array2table(Bt_range,'RowNames',rowNames,'VariableNames',colNames)
     
    writetable(pc_change_table,'pc_change.csv')
     
    % Now plot a separate chart for each of the coefficients
     
    figure
    for i = 1:21       
        subplot(7,3,i)
        plot(1:t,Bt_postmean(i,:))  
    end 
     
    % Now we move to plotting the IRF. This section takes moments along the
    % first dimension, i.e. across the Gibbs sample iterations. The moments
    % are for the 16th, 50th and 84th percentile. 
     
        qus = [.16, .5, .84];
        imp75XY=squeeze(quantile(imp75,qus)); 
        imp81XY=squeeze(quantile(imp81,qus));
        imp96XY=squeeze(quantile(imp96,qus));
        
        % Plot impulse responses
        figure       
        set(0,'DefaultAxesColorOrder',[0 0 0],...
            'DefaultAxesLineStyleOrder','--|-|--')
        subplot(3,3,1)
        plot(1:nhor,squeeze(imp75XY(:,1,:)))
        title('Impulse response of inflation, 1975:Q1')
        xlim([1 nhor])
        ylim([-0.2 0.1])
    %    % yline(0)
        set(gca,'XTick',0:3:nhor)
        subplot(3,3,2)
        plot(1:nhor,squeeze(imp75XY(:,2,:)))
        title('Impulse response of unemployment, 1975:Q1')
        xlim([1 nhor])
        ylim([-0.2 0.2])
       % yline(0)
        set(gca,'XTick',0:3:nhor)    
        subplot(3,3,3)
        %ylim([0 1])
    %    % yline(0)
        plot(1:nhor,squeeze(imp75XY(:,3,:)))
        title('Impulse response of interest rate, 1975:Q1')
        xlim([1 nhor])
        %ylim([-0.3 0.1])
       % yline(0)
        set(gca,'XTick',0:3:nhor)    
        subplot(3,3,4)
        plot(1:nhor,squeeze(imp81XY(:,1,:)))
        title('Impulse response of inflation, 1981:Q3')
        xlim([1 nhor])
        ylim([-0.2 0.1])
       % yline(0)
        set(gca,'XTick',0:3:nhor)    
        subplot(3,3,5)
        plot(1:nhor,squeeze(imp81XY(:,2,:)))
        title('Impulse response of unemployment, 1981:Q3')
        xlim([1 nhor])
        ylim([-0.2 0.2])
       % yline(0)
        set(gca,'XTick',0:3:nhor)    
        subplot(3,3,6)
        plot(1:nhor,squeeze(imp81XY(:,3,:)))
        title('Impulse response of interest rate, 1981:Q3')
        xlim([1 nhor])
        %ylim([-0.4 0.1])
       % yline(0)
        set(gca,'XTick',0:3:nhor)    
        subplot(3,3,7)
        plot(1:nhor,squeeze(imp96XY(:,1,:)))
        title('Impulse response of inflation, 1996:Q1')
        xlim([1 nhor])
        ylim([-0.2 0.1])
       % yline(0)
        set(gca,'XTick',0:3:nhor)    
        subplot(3,3,8)
        plot(1:nhor,squeeze(imp96XY(:,2,:)))
        title('Impulse response of unemployment, 1996:Q1')
        xlim([1 nhor])
        ylim([-0.2 0.2])
       % yline(0)
        set(gca,'XTick',0:3:nhor)
        subplot(3,3,9)
        plot(1:nhor,squeeze(imp96XY(:,3,:)))
        title('Impulse response of interest rate, 1996:Q1')
        xlim([1 nhor])
         %ylim([0 1])
        % yline(0)
        set(gca,'XTick',0:3:nhor)
        
    disp('             ')
    disp('To plot impulse responses, use:         plot(1:nhor,squeeze(imp75XY(:,VAR,:)))           ')
    disp('             ')
    disp('where VAR=1 for impulses of inflation, VAR=2 for unemployment and VAR=3 for interest rate')
    

    三、运行结果

    在这里插入图片描述
    在这里插入图片描述

    四、备注

    版本:2014a

    展开全文
  • matlab自相关代码资料夹结构 Data_stock_market_return.mat 它包含最终数据集。 股市指数的观察来自https://fr.finance.yahoo.com/ 。 该研究的数据集包含每个系列的3514个有效数据点(从2002年12月16日到2017年12月...
  • matlab 在险价值 VaR 的计算

    千次阅读 多人点赞 2019-04-17 12:13:26
    matlab 在险价值 VaR 的计算 matlab 在险价值 VaR 的计算 VaR 模型 数据获取 历史模拟法 蒙特卡罗模拟法 参数模型法 代码和数据下载 VaR 模型 Value at Risk 在险价值,即 VaR。是指一定时期内,一定置信...

    matlab 在险价值 VaR 的计算

    VaR 模型

    Value at Risk 在险价值,即 VaR。是指一定时期内,一定置信水平下,某种资产组合面临的最大损失。

    Prob(ΔpVaR)=1αProb(\Delta p \le VaR) = 1 - \alpha

    Δp\Delta p 是该时期内投资组合市值变动。这个式子表示,在一定持有期 Δt\Delta t 内,在给定的置信水平 1α1- \alpha 下,该投资组合的最大损失不会超过 VaR。

    数据获取

    以前推荐用tushare是因为tushare免费,但是现在tushare非常无耻,积分不透明变像收费早已背离初心,建议用其他数据平台。

    在做 VaR 实战之前,我们得先获取实验数据。这部分可以参考我之前的文章。简单来说,就是需要个 tuhsare 账号,可以点去这里注册个https://tushare.pro/register?reg=126259账号,好像直接点击链接会提示风险,需要复制打地址栏打开,然后获取 自己的token,替换下文中的 token。代码如下:

    close all
    clear
    clc
    %%  数据准备
    % 尝试用从tushare读取数据,如果因为版本问题读取数据失败
    % 则用之前保存下来的数据进行计算
    try
        % 加载tushare包
        addpath(genpath(pwd));
        % 替换成你自己的token
        token = '**************937459c8b611e0a97d0abf*******';
        api = pro_api(token);
        stockall = api.query('stock_basic');
        % 读取数据参数设定
        start_time = '20180101';
        end_time = '';
        ktype = 'D';
        % 取沪深300做为市场指数
        indexdata = pro_bar('000300.SH', api, start_time, end_time,ktype,'I');
        indexdata = flipud(indexdata);
        % 取3支股票
        nstock = 3;
        % 记录用到的3支股票
        stocklist = [];
        stockdata = cell(nstock,1);
        nday = size(indexdata,1);
        % closeprice 第一列为指数价格,其他列为股票数据
        closeprice = indexdata.close;
        temp = indexdata.trade_date;
        temp = char(temp);
        temp = str2num(temp);
        tradedate = datetime(temp,'ConvertFrom','yyyymmdd','format','yyyy-MM-dd');
        
        % 获取股票数据
        m = 0;
        for i = 1:size(stockall,1)
            temp = pro_bar(stockall.ts_code{i}, api, start_time, end_time,ktype,'E','qfq');
            % 有的股票有停牌,我们得选取没有停牌的股票
            if size(temp,1) == nday
                m = m+1;
                stockdata{m} = flipud(temp);
                closeprice(:,m+1) = stockdata{m}.close;
                stocklist = [stocklist;stockall(i,:)];
            end
            if size(stocklist,1) == nstock
                break;
            end
        end
        save('tempdata.mat');
    catch
        load tempdata
    end
    
    % 对所选用的股票进行绘图
    plot(tradedate,closeprice(:,2))
    hold on
    plot(tradedate,closeprice(:,3))
    plot(tradedate,closeprice(:,4))
    xlabel('时间');
    ylabel('股价');
    legend(stocklist.name);
    title('投资组合成分股股价变动图');
    saveas(gcf,'投资组合成分股股价变动图.jpg');
    

    这部分的代码我进行了容错处理,如果因为版本之类的原因,你的 matlb 无法正确运行这部分的代码,可以去文末我附的链接中下载,这部分的代码和数据,这样也能正确运行代码。这部分绘制出图像:

    投资组合成分股股价变动图.jpg

    历史模拟法

    所以历史模拟法,就是用投资组合的历史收益率经验分布的分位数来计算 VaR。代码如下:

    %% 历史模拟法
    %  假定每只股票买N
    N = 1000;
    % 假定投资组合等权重,求投资组合市值
    value = (closeprice(:,2) + closeprice(:,3) + closeprice(:,4)) * N;
    % 投资组合收益率
    ret = price2ret(value);
    figure;
    subplot(2,1,1);
    plot(tradedate,value);
    xlabel('时间');
    ylabel('组合市值');
    title('投资组合市值')
    subplot(2,1,2);
    plot(tradedate(2:end),ret,'*');
    xlabel('时间');
    ylabel('收益率');
    title('投资组合日收益率')
    saveas(gcf,'投资组合市值及收益率.jpg');
    %绘制投资组合收益率直方图
    figure;
    histogram(ret,20);
    ylabel('天数');
    xlabel('投资组合日收益率');
    title('历史模拟法投资组合日收益率直方图');
    saveas(gcf,'历史模拟法投资组合日收益率直方图.jpg')
    %在5%置信度时,市值亏损的最大比率
    Var = -prctile(ret,5) * value(end);
    disp(['历史模拟法投资组合VaR为',num2str(Var)]);
    
    

    首先,这边得解释一下投资组合的构建。这里是假定,每只股票持有 1000 股。然后用 prctile 函数求出经验分布的分位数。这部分结果如下:

    投资组合市值及收益率.jpg

    历史模拟法投资组合日收益率直方图.jpg

    命令行中输出:历史模拟法投资组合VaR为1860.6194

    蒙特卡罗模拟法

    蒙特卡罗法也就是根据历史的收益率估计出均值和波动率,然后再用这个均值和波动率,模拟出未来价格路径。在分别求每条价格路径的 VaR ,最后求其均值。若是对蒙特卡罗法本身有疑问,可以简单参考一下这篇文章

    %% 蒙特卡罗模拟
    %计算日收益率均值和方差
    mu = mean(ret);
    vol = std(ret);
    % 年化波动率
    vol = vol * sqrt(250);
    %资产初始值
    s0 = value(end);
    % 模拟时长
    T = 1;
    % 模拟间隔点
    nStep = 250;
    % 模拟路径数
    nPath = 1000;
    %用蒙特卡洛模拟1000次
    sPath = simulatePath(s0,mu,vol,T,nStep,nPath);
    % 绘制蒙特卡罗法模拟路径
    figure;
    plot(sPath);
    xlabel('模拟时间点');
    ylabel('组合市值');
    title('蒙特卡罗法模拟路径');
    saveas(gcf,'蒙特卡罗法模拟路径.jpg')
    % 计算每条模拟路径的var
    ret_mc = price2ret(sPath);
    Var = -mean(prctile(ret_mc,5)) * value(end);
    disp(['蒙特卡罗法投资组合VaR为',num2str(Var)]);
    

    其中 simulatePath 函数如下:

    function sPath = simulatePath(S0,mu,sigma,T,nStep,nPath)
    % mu=r-q
    
    % time step
    deltaT = T / nStep;
    % initialize stock price matrix
    sPath = nan(nStep+1,nPath);
    sPath(1,:) = S0;
    % Simulated according to ITO's lemma
    p1 = (mu - 0.5 * sigma ^ 2 ) * deltaT;
    p2 = sigma * sqrt(deltaT);
    for i = 1:nPath
        for j = 1:nStep
            sPath(j+1,i) = sPath(j,i) * exp(p1 + p2 * randn);
        end
    end
    
    
    end
    

    蒙特卡罗法模拟路径.jpg

    命令行输出:蒙特卡罗法投资组合VaR为1816.0601

    参数模型法

    所谓参数模型法就是调用了 matlab 内置函数 portvrisk。

    %% 参数模型法
    vol = std(ret);
    ValueAtRisk = portvrisk(mu,vol);
    Var = ValueAtRisk * value(end);
    disp(['参数模型法投资组合VaR为',num2str(Var)]);
    

    命令行输出:参数模型法投资组合VaR为1811.7731

    代码和数据下载

    我已经将本文的代码和数据上传至 csdn资源下载页面,欢迎下载。随便骂一句,垃圾csdn代码高亮都搞不清楚,老子迟早用自己博客。

    展开全文
  • 此例程将向量自回归 (VAR) 的参数估计映射到相应移动平均 (MA) 模型的参数估计中。 此函数的输出可用于构建 VAR 模型的结构脉冲响应函数。
  • 向量的散度matlab代码短信增值业务 用于航空系统异常检测的半马尔可夫切换矢量自回归 (SMS-VAR) 模型 用于生成、训练和评估 SMS-VAR 模型的代码演示版。 运行 Matlab 脚本 runSim.m 来执行代码。 (确保包含代码的...
  • matlab分时代码BVAR连接 描述 这是一种用户友好的Matlab GUI,它对贝叶斯多主题向量自回归(VAR模型实施了变分推理方法,以便基于静止状态功能MRI数据来推理有效的大脑连通性。 建模框架使用贝叶斯变量选择方法,...
  • 广义泊松回归代码matlab实现
  • 一、 VAR:向量自回归模型,结果仅具有统计上的 意义 SVAR:结构向量自回归模型 TVP-VAR:Time Varying Parameter-Stochastic Volatility-Vector Auto Regression。...在VAR模型中,所有的参数遵循一阶游走过程

    一、简介

    VAR:向量自回归模型,结果仅具有统计上的 意义

    SVAR:结构向量自回归模型

    TVP-VAR:Time Varying Parameter-Stochastic Volatility-Vector Auto Regression。时变参数随机波动率向量自回归模型,与VAR 不同的是,模型没有同方差的假定,更符合实际。并且时变参数假定随机波动率,更能捕捉到经济变量在不同时代背景下所具有的关系和特征(时变影响)。将随机波动性纳入TVP估算中可以显着提高估算性能。在VAR模型中,所有的参数遵循一阶游走过程。

      随机波动的概念在TVP-VAR中很重要,随机波动是1976年由Black提出,随后在经济计量中有很大的发展。近几年,随机波动也经常被用在宏观经济的经验分析中。很多情况下,经济数据的产生过程中具有漂移系数和随机波动的冲击,如果是这种情况,那么使用具有时变系数但具有恒定波动性的模型会引起一个问题,即由于忽略了扰动中波动性的可能变化,估计的时变系数可能会出现偏差。为了避免这种问题,TVP-VAR模型假设了随机波动性。尽管似然函数难以处理,所以随机波动率使估计变得困难,但是可以在贝叶斯推断中使用马尔可夫链蒙特卡罗(MCMC)方法来估计模型。
    

    从 具有随机波动性的时变参数(TVP)回归模型的估计算法(TVP-VAR模型的单变量情况)说起。
    在这里插入图片描述
    在这里插入图片描述
    二、MCMC算法

    MCMC算法,一种随机采样方法,即蒙特卡罗方法(Monte Carlo Simulation,简称MC)和马尔科夫链(Markov Chain ,也简称MC)。MCMC方法是用来在概率空间,通过随机采样估算兴趣参数的后验分布。

    MCMC方法是在贝叶斯推断的背景下考虑的,其目标是在研究人员预先设定的特定先验概率密度下评估目标参数的联合后验分布。

    在贝叶斯推断中,为未知参数θ的向量指定先验密度,用π(θ)表示。用f(y|θ)表示数据y = {y 1 ,…,y n }的似然函数。然后根据先验分布进行推断,用 π(θ|y)表示。

    在这里插入图片描述
    MCMC算法通过递归采样条件后验分布来进行,其中在仿真中使用了条件参数的最新值。

    Gibbs采样器是著名的MCMC方法之一。步骤如下:
    在这里插入图片描述
    对于MCMC算法,π(θ)为先验密度,后验分布 π(θ,α,h|y)^4。
    在这里插入图片描述
    (具体参数参考Nakajima原文)

    三、TVP-VAR Model

    首先考虑SVAR model的定义:
    在这里插入图片描述
    yt是k×1的观察变量,A,F是k×k的系数矩阵,扰动项ut是k×1的结构性冲击,假设,
    在这里插入图片描述
    假设A为下三角矩阵,通过递归识别指定结构冲击的关系。

    在这里插入图片描述

    所以,式6可以写为:
    在这里插入图片描述
    在这里插入图片描述
    3.当在贝叶斯推断中实现TVP-VAR模型时,应谨慎选择先验,因为TVP-VAR模型具有许多状态变量,并且其过程被建模为非平稳随机游走过程TVP-VAR模型非常灵活,状态变量可以捕捉潜在经济结构的渐进和突变。但是在VAR模型中的每个参数中允许时间变化可能会导致过度识别问题。

    四、进行TVP-VAR建模时,也需要数据平稳。可以用ADF单位根检验法检验数据的平稳性,不平稳的数据可以做差分,直至平稳,用差分后的数据进行建模。

    需要注意的是:TVP-VAR建模时,变量之间的顺序会影响到最后的实证结果。变量的顺序问题可能是实证结果不符合预期的一个原因。最好把关注的变量放在首位,这样在时变关系图中就能得到较好的展示。

    二、源代码

    %==================================4.1 Main================================
    clear;
    clc;
     
    % Load Korobilis (2008) quarterly data
    load ydata.dat; % data 
    load yearlab.dat; % data labels 
     
    %%
    %----------------------------------BASICS----------------------------------
     
     
    Y=ydata;
    t=size(Y,1);        % t - The total number of periods in the raw data (t=215)
    M=size(Y,2);        % M - The dimensionality of Y (i.e. the number of variables)(M=3)
    tau = 40;           % tau - the size of the training sample (the first forty quarters)
    p = 2;              % p - number of lags in the VAR model 
     
    %% Generate the Z_t matrix, i.e. the regressors in the model. 
     
    ylag = mlag2(Y,p);          % This function generates a 215x6 matrix with p lags of variable Y. 
    ylag = ylag(p+tau+1:t,:);   % Then remove our training sample, so now a 173x6 matrix. 
     
    K = M + p*(M^2);            % K is the number of elements in the state vector
     
    % Here we distribute the lagged y data into the Z matrix so it is
    % conformable with a beta_t matrix of coefficients. 
     
    Z = zeros((t-tau-p)*M,K);   
    for i = 1:t-tau-p
        ztemp = eye(M);
        for j = 1:p        
            xtemp = ylag(i,(j-1)*M+1:j*M); 
            xtemp = kron(eye(M),xtemp);
            ztemp = [ztemp xtemp];
        end
        Z((i-1)*M+1:i*M,:) = ztemp;
    end
     
    % Redefine our variables to exclude the training sample and the first two
    % lags that we take as given, taking total number of periods (t) from 215 
    % to 173. 
     
    y = Y(tau+p+1:t,:)';
    yearlab = yearlab(tau+p+1:t);
    t=size(y,2); % t now equals 173
     
    %% --------------------MODEL AND GIBBS PRELIMINARIES-----------------------
     
    nrep = 5000;    % Number of sample draws 
    nburn = 2000;   % Number of burn-in-draws
    it_print = 100; % Print in the screen every "it_print"-th iteration
     
     
    %% INITIAL STATE VECTOR PRIOR  
     
    % We use the first 40 observations (tau) to run a standard OLS of the
    % measurement equation, using the function ts_prior. The result is
    % estimates for priors for B_0 and Var(B_0). 
     
    [B_OLS,VB_OLS]= ts_prior(Y,tau,M,p);
     
    % Given the distributions we have, we now have to define our priors for B,
    % Q and Sigma. These are set in accordance with how they are set in 
    % Primiceri (2005). These are the hyperparameters of the beta, Q and Sigma
    % initial priors. 
     
    B_0_prmean = B_OLS;
    B_0_prvar = 4*VB_OLS; 
     
    Q_prmean = ((0.01).^2)*tau*VB_OLS;
    Q_prvar = tau;
     
    Sigma_prmean = eye(M);
    Sigma_prvar = M+1;
     
    % To start the Kalman filtering assign arbitrary values that are in support
    % of their priors, Q and Sigma. 
     
    consQ = 0.0001;
    Qdraw = consQ*eye(K);
    Sigmadraw = 0.1*eye(M);
     
    % Create some matrices for storage that will be filled in once we
    % start the Gibbs sampling.
     
    Btdraw = zeros(K,t); 
    Bt_postmean = zeros(K,t);
    Qmean = zeros(K,K);
    Sigmamean = zeros(M,M);
     
    %% -------------------------IRF-PRELIMINARIES------------------------------
    nhor = 21;     % The number of periods in the impulse response function. 
     
    % Matricies to be filled containing IRFs for 1975q1, 1981q3, 1996q1. The 
    % dimensions correspond to the iterations of the gibbs sample, each of the
    % variables, and each of the 21 periods of the IRF analysis. 
     
    imp75 = zeros(nrep,M,nhor); 
    imp81 = zeros(nrep,M,nhor);
    imp96 = zeros(nrep,M,nhor);
     
    % This corresponds to variable J introduced in equation (14) in the report
     
    bigj = zeros(M,M*p); 
    bigj(1:M,1:M) = eye(M);
     
     
    %% ================ START GIBBS SAMPLING ==================================
     
    tic; % This is just a timer
    disp('Number of iterations');
     
     
    for irep = 1:nrep + nburn    % 7000 gibbs iterations starts here
        % Print iterations - this just updates on the progress of the sampling
        if mod(irep,it_print) == 0
            disp(irep);toc;
        end
        
    %% Draw 1: B_t from p(B_t|y,Sigma)
        
        % We use the function 'carter_kohn_hom' to to run the FFBS algorithm. 
        % This results in a 21x173 matrix, corresponding to one Gibbs sample
        % draw of each of the coefficients in each time period. The inputs
        % Sigmadraw and Qdraw are updated for each Gibbs sample repetition.
        
        [Btdraw] = carter_kohn_hom(y,Z,Sigmadraw,Qdraw,K,M,t,B_0_prmean,B_0_prvar);
        
      
    %% Draw 2: Q from p(Q^{-1}|y,B_t) which is i-Wishart
     
        % We draw Q from an Inverse Wishart distribution. The parameters 
        % of the distribution are derived as equation (11) in the main report.
        % The mean is taken as the inverse of the accumulated sum of squared 
        % errors added to the prior mean, and the variance is simply t.  
        
        % Differencing Btdraw to create the sum of squared errors
        
        Btemp = Btdraw(:,2:t)' - Btdraw(:,1:t-1)'; 
        sse_2Q = zeros(K,K);
        for i = 1:t-1
            sse_2Q = sse_2Q + Btemp(i,:)'*Btemp(i,:);
        end
     
        Qinv = inv(sse_2Q + Q_prmean);      % compute mean to use for Wishart draw
        Qinvdraw = wish(Qinv,t+Q_prvar);    % draw inv q from the wishart distribution
        Qdraw = inv(Qinvdraw);              % find non-inverse q 
        
    %% Draw 3: Sigma from p(Sigma|y,B_t) which is i-Wishart
     
        % We draw Sigma from an Inverse Wishart distribution. The parameters 
        % of the distirbution are derived as equation (10) in the main report. 
        % The mean is taken as the inverse of the sum of squared residuals
        % added to the prior mean. The variance is simply t. 
        
        % Find residuals using data and the current draw of coefficients
        
        resids = zeros(M,t);
        for i = 1:t
            resids(:,i) = y(:,i) - Z((i-1)*M+1:i*M,:)*Btdraw(:,i);
        end
        
        % Create a matrix for the accumulated sum of squared residuals, to
        % be used as the mean parameter in the i-Wishart draw below. 
        sse_2S = zeros(M,M);
        for i = 1:t
            sse_2S = sse_2S + resids(:,i)*resids(:,i)';
        end
        
        Sigmainv = inv(sse_2S + Sigma_prmean);          % compute mean to use for the Wishart
        Sigmainvdraw = wish(Sigmainv,t+Sigma_prvar);    % draw from the Wishsart distribution
        Sigmadraw = inv(Sigmainvdraw);                  % turn into non-inverse Sigma
        
    %% IRF 
        % We only apply IRF analysis once we have exceeded the burn-in draws phase.
        if irep > nburn;         
                
                % Create matrix that is going to contain all beta draws over
                % which we will take the mean after the Gibbs sampler as our moment
                % estimate: 
                Bt_postmean = Bt_postmean + Btdraw;
                
                % biga is the A matrix of the VAR(1) version of our VAR(2) model, 
                % found in equation (12. biga changes in every period of the 
                % analysis, because the coefficients are time varying, so we 
                % apply the analysis below in every time period. 
                
                biga = zeros(M*p,M*p); 
                for j = 1:p-1
                    biga(j*M+1:M*(j+1),M*(j-1)+1:j*M) = eye(M); % fill the A matrix with identity matrix (3) in bottom left corner
                end
     
                % The following procedure is applied separately in each time period. 
                
                % This loop takes coefficients of the relevant time period from
                % Bt_draw (which contains all coefficients for all t) and uses
                % them to update the biga matrix, so that it can change for
                % every t. 
                
                for i = 1:t 
                    bbtemp = Btdraw(M+1:K,i);  % get the draw of B(t) at time i=1,...,T  (exclude intercept)
                    splace = 0;
                    for ii = 1:p
                        for iii = 1:M
                            biga(iii,(ii-1)*M+1:ii*M) = bbtemp(splace+1:splace+M,1)'; % Load non-intercept coefficient draws
                            splace = splace + M;
                        end
                    end
                    
                    % Next we want to create a shock matrix in which the third
                    % column is [0 0 1]', therefore implementing a unit shock
                    % in the interest rate. 
                    
                    shock = eye(3);
                    
                    % Now get impulse responses for 1 through nhor future
                    % periods. impresp is a 3x63 matrix which contains 9
                    % response values in total for each period, 3 for each 
                    % variable. These three responses correspond to the 3
                    % possible shocks that are contained in the schock
                    % matrix
                    
                    % bigai is updated through mulitiplication with the 
                    % coefficient matrix after each time period. 
                                    
                    % Create a results matrix to store impulse responses in all periods
                    impresp = zeros(M,M*nhor); 
                    
                    % Fill in the first period of the results matrix with the shock (as defined above) 
                    impresp(1:M,1:M) = shock;
                    
                    % Create a separate variable for the a matrix so that we
                    % can update it for each period of the IRF analysis. 
                    bigai = biga; 
                    
                    % This follows the impulse response function as in equation 15.
                    % Fill in each period of the results matrix according to
                    % the impulse response function formula. 
                    
                    for j = 1:nhor-1
                        impresp(:,j*M+1:(j+1)*M) = bigj*bigai*bigj'*shock;
                        bigai = bigai*biga; % update the coefficient matrix for next period
                    end
     
                    % The section below keeps only the responses that we are interested in:
                    % - those from the periods 1975q1, 1981q3, and 1996q1
                    % - those that correspond to the shock in the interest 
                    % rate (i.e. those caused by the third column of our shock
                    % matrix). 
                    
                    if yearlab(i,1) == 1975.00;   % store only IRF from 1975:Q1
                        impf_m = zeros(M,nhor);
                        jj=0;
                        for ij = 1:nhor
                            jj = jj + M;    % select only the third column for each time period of the IRF
                            impf_m(:,ij) = impresp(:,jj);
                        end           
                    % For each iteration of the Gibbs sampler, fill in the
                    % results along the first dimension 
                        imp75(irep-nburn,:,:) = impf_m; 
                    end
                    
                    if yearlab(i,1) == 1981.50;   % store only IRF from 1981:Q3
                        impf_m = zeros(M,nhor);
                        jj=0;
                        for ij = 1:nhor
                            jj = jj + M;    % select only the third column for each time period of the IRF
                            impf_m(:,ij) = impresp(:,jj);
                        end
                    % For each iteration of the Gibbs sample, fill in the
                    % results along the first dimension 
                        imp81(irep-nburn,:,:) = impf_m; 
                    end
                    
                    if yearlab(i,1) == 1996.00;   % store only IRF from 1996:Q1
                        impf_m = zeros(M,nhor);
                        jj=0;
                        for ij = 1:nhor
                            jj = jj + M;    % select only the third column for each time period of the IRF
                            impf_m(:,ij) = impresp(:,jj);
                        end
                    % For each iteration of the Gibbs sample, fill in the
                    % results along the first dimension 
                        imp96(irep-nburn,:,:) = impf_m; 
                    end
                end % End getting impulses for each time period 
            end % End the impulse response calculation section   
    end % End main Gibbs loop (for irep = 1:nrep+nburn)
     
    clc;
    toc; % Stop timer and print total time
    %% ================ END GIBBS SAMPLING ==================================
     
    % Even though it is not used in our IRF analysis since we are integrating
    % that into the Gibbs sampling loop, here is how to take the mean of the 
    % draw of the betas as moment estimate:
     
    Bt_postmean = Bt_postmean./nrep;
     
    %% Graphs and tables
     
    % This works out the percentage range of that each variable's coefficient
    % spans across time. I.e. to what extent was the TVP facility used by each
    % variable in the model? This is calculated by finding the range for each
    % variable as a percentage of the mean coefficient size for that variable.
    % The result is a 21x1 vector, and it is reshaped into a 3x7 matrix in order
    % to map onto the system of equations (2,3,and 4) set out in the report.
     
    Bt_range = ones(21,1)
    for i = 1:21
        Bt_range(i) = abs((max(Bt_postmean(i,:))-min(Bt_postmean(i,:)))/mean(Bt_postmean(i,:)))
    end 
    Bt_range = reshape(Bt_range,3,7)
     
    % Create a table of coefficient ranges for export to the report
     
    rowNames = {'Inflation','Unemployment','Interest Rate'};
    colNames = {'Intercept','Inf_1','Unemp_1', 'IR_1','Inf_2','Unemp_2', 'IR_2'};
    pc_change_table = array2table(Bt_range,'RowNames',rowNames,'VariableNames',colNames)
     
    writetable(pc_change_table,'pc_change.csv')
     
    % Now plot a separate chart for each of the coefficients
     
    figure
    for i = 1:21       
        subplot(7,3,i)
        plot(1:t,Bt_postmean(i,:))  
    end 
     
    % Now we move to plotting the IRF. This section takes moments along the
    % first dimension, i.e. across the Gibbs sample iterations. The moments
    % are for the 16th, 50th and 84th percentile. 
     
        qus = [.16, .5, .84];
        imp75XY=squeeze(quantile(imp75,qus)); 
        imp81XY=squeeze(quantile(imp81,qus));
        imp96XY=squeeze(quantile(imp96,qus));
        
        % Plot impulse responses
        figure       
        set(0,'DefaultAxesColorOrder',[0 0 0],...
            'DefaultAxesLineStyleOrder','--|-|--')
        subplot(3,3,1)
        plot(1:nhor,squeeze(imp75XY(:,1,:)))
        title('Impulse response of inflation, 1975:Q1')
        xlim([1 nhor])
        ylim([-0.2 0.1])
    %    % yline(0)
        set(gca,'XTick',0:3:nhor)
        subplot(3,3,2)
        plot(1:nhor,squeeze(imp75XY(:,2,:)))
        title('Impulse response of unemployment, 1975:Q1')
        xlim([1 nhor])
        ylim([-0.2 0.2])
       % yline(0)
        set(gca,'XTick',0:3:nhor)    
        subplot(3,3,3)
        %ylim([0 1])
    %    % yline(0)
        plot(1:nhor,squeeze(imp75XY(:,3,:)))
        title('Impulse response of interest rate, 1975:Q1')
        xlim([1 nhor])
        %ylim([-0.3 0.1])
       % yline(0)
        set(gca,'XTick',0:3:nhor)    
        subplot(3,3,4)
        plot(1:nhor,squeeze(imp81XY(:,1,:)))
        title('Impulse response of inflation, 1981:Q3')
        xlim([1 nhor])
        ylim([-0.2 0.1])
       % yline(0)
        set(gca,'XTick',0:3:nhor)    
        subplot(3,3,5)
        plot(1:nhor,squeeze(imp81XY(:,2,:)))
        title('Impulse response of unemployment, 1981:Q3')
        xlim([1 nhor])
        ylim([-0.2 0.2])
       % yline(0)
        set(gca,'XTick',0:3:nhor)    
        subplot(3,3,6)
        plot(1:nhor,squeeze(imp81XY(:,3,:)))
        title('Impulse response of interest rate, 1981:Q3')
        xlim([1 nhor])
        %ylim([-0.4 0.1])
       % yline(0)
        set(gca,'XTick',0:3:nhor)    
        subplot(3,3,7)
        plot(1:nhor,squeeze(imp96XY(:,1,:)))
        title('Impulse response of inflation, 1996:Q1')
        xlim([1 nhor])
        ylim([-0.2 0.1])
       % yline(0)
        set(gca,'XTick',0:3:nhor)    
        subplot(3,3,8)
        plot(1:nhor,squeeze(imp96XY(:,2,:)))
        title('Impulse response of unemployment, 1996:Q1')
        xlim([1 nhor])
        ylim([-0.2 0.2])
       % yline(0)
        set(gca,'XTick',0:3:nhor)
        subplot(3,3,9)
        plot(1:nhor,squeeze(imp96XY(:,3,:)))
        title('Impulse response of interest rate, 1996:Q1')
        xlim([1 nhor])
         %ylim([0 1])
        % yline(0)
        set(gca,'XTick',0:3:nhor)
        
    disp('             ')
    disp('To plot impulse responses, use:         plot(1:nhor,squeeze(imp75XY(:,VAR,:)))           ')
    disp('             ')
    disp('where VAR=1 for impulses of inflation, VAR=2 for unemployment and VAR=3 for interest rate')
    

    三、运行结果

    在这里插入图片描述
    在这里插入图片描述

    四、备注

    版本:2014a

    展开全文
  • 本资源包含,用matlab实现历史模拟法、蒙特卡罗法、参数模型法等三种方法求解VaR
  • 一、 VAR:向量自回归模型,结果仅具有统计上的 意义 SVAR:结构向量自回归模型 TVP-VAR:Time Varying Parameter-Stochastic Volatility-Vector Auto Regression。...在VAR模型中,所有的参数遵循一阶游走过程

    一、简介

    VAR:向量自回归模型,结果仅具有统计上的 意义

    SVAR:结构向量自回归模型

    TVP-VAR:Time Varying Parameter-Stochastic Volatility-Vector Auto Regression。时变参数随机波动率向量自回归模型,与VAR 不同的是,模型没有同方差的假定,更符合实际。并且时变参数假定随机波动率,更能捕捉到经济变量在不同时代背景下所具有的关系和特征(时变影响)。将随机波动性纳入TVP估算中可以显着提高估算性能。在VAR模型中,所有的参数遵循一阶游走过程。

      随机波动的概念在TVP-VAR中很重要,随机波动是1976年由Black提出,随后在经济计量中有很大的发展。近几年,随机波动也经常被用在宏观经济的经验分析中。很多情况下,经济数据的产生过程中具有漂移系数和随机波动的冲击,如果是这种情况,那么使用具有时变系数但具有恒定波动性的模型会引起一个问题,即由于忽略了扰动中波动性的可能变化,估计的时变系数可能会出现偏差。为了避免这种问题,TVP-VAR模型假设了随机波动性。尽管似然函数难以处理,所以随机波动率使估计变得困难,但是可以在贝叶斯推断中使用马尔可夫链蒙特卡罗(MCMC)方法来估计模型。
    

    从 具有随机波动性的时变参数(TVP)回归模型的估计算法(TVP-VAR模型的单变量情况)说起。
    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述
    对于MCMC算法,π(θ)为先验密度,后验分布 π(θ,α,h|y)^4。
    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述
    在这里插入图片描述
    四、进行TVP-VAR建模时,也需要数据平稳。可以用ADF单位根检验法检验数据的平稳性,不平稳的数据可以做差分,直至平稳,用差分后的数据进行建模。

    需要注意的是:TVP-VAR建模时,变量之间的顺序会影响到最后的实证结果。变量的顺序问题可能是实证结果不符合预期的一个原因。最好把关注的变量放在首位,这样在时变关系图中就能得到较好的展示。

    二、源代码

    %==================================4.1 Main================================
    clear;
    clc;
     
    % Load Korobilis (2008) quarterly data
    load ydata.dat; % data 
    load yearlab.dat; % data labels 
     
    %%
    %----------------------------------BASICS----------------------------------
     
     
    Y=ydata;
    t=size(Y,1);        % t - The total number of periods in the raw data (t=215)
    M=size(Y,2);        % M - The dimensionality of Y (i.e. the number of variables)(M=3)
    tau = 40;           % tau - the size of the training sample (the first forty quarters)
    p = 2;              % p - number of lags in the VAR model 
     
    %% Generate the Z_t matrix, i.e. the regressors in the model. 
     
    ylag = mlag2(Y,p);          % This function generates a 215x6 matrix with p lags of variable Y. 
    ylag = ylag(p+tau+1:t,:);   % Then remove our training sample, so now a 173x6 matrix. 
     
    K = M + p*(M^2);            % K is the number of elements in the state vector
     
    % Here we distribute the lagged y data into the Z matrix so it is
    % conformable with a beta_t matrix of coefficients. 
     
    Z = zeros((t-tau-p)*M,K);   
    for i = 1:t-tau-p
        ztemp = eye(M);
        for j = 1:p        
            xtemp = ylag(i,(j-1)*M+1:j*M); 
            xtemp = kron(eye(M),xtemp);
            ztemp = [ztemp xtemp];
        end
        Z((i-1)*M+1:i*M,:) = ztemp;
    end
     
    % Redefine our variables to exclude the training sample and the first two
    % lags that we take as given, taking total number of periods (t) from 215 
    % to 173. 
     
    y = Y(tau+p+1:t,:)';
    yearlab = yearlab(tau+p+1:t);
    t=size(y,2); % t now equals 173
     
    %% --------------------MODEL AND GIBBS PRELIMINARIES-----------------------
     
    nrep = 5000;    % Number of sample draws 
    nburn = 2000;   % Number of burn-in-draws
    it_print = 100; % Print in the screen every "it_print"-th iteration
     
     
    %% INITIAL STATE VECTOR PRIOR  
     
    % We use the first 40 observations (tau) to run a standard OLS of the
    % measurement equation, using the function ts_prior. The result is
    % estimates for priors for B_0 and Var(B_0). 
     
    [B_OLS,VB_OLS]= ts_prior(Y,tau,M,p);
     
    % Given the distributions we have, we now have to define our priors for B,
    % Q and Sigma. These are set in accordance with how they are set in 
    % Primiceri (2005). These are the hyperparameters of the beta, Q and Sigma
    % initial priors. 
     
    B_0_prmean = B_OLS;
    B_0_prvar = 4*VB_OLS; 
     
    Q_prmean = ((0.01).^2)*tau*VB_OLS;
    Q_prvar = tau;
     
    Sigma_prmean = eye(M);
    Sigma_prvar = M+1;
     
    % To start the Kalman filtering assign arbitrary values that are in support
    % of their priors, Q and Sigma. 
     
    consQ = 0.0001;
    Qdraw = consQ*eye(K);
    Sigmadraw = 0.1*eye(M);
     
    % Create some matrices for storage that will be filled in once we
    % start the Gibbs sampling.
     
    Btdraw = zeros(K,t); 
    Bt_postmean = zeros(K,t);
    Qmean = zeros(K,K);
    Sigmamean = zeros(M,M);
     
    %% -------------------------IRF-PRELIMINARIES------------------------------
    nhor = 21;     % The number of periods in the impulse response function. 
     
    % Matricies to be filled containing IRFs for 1975q1, 1981q3, 1996q1. The 
    % dimensions correspond to the iterations of the gibbs sample, each of the
    % variables, and each of the 21 periods of the IRF analysis. 
     
    imp75 = zeros(nrep,M,nhor); 
    imp81 = zeros(nrep,M,nhor);
    imp96 = zeros(nrep,M,nhor);
     
    % This corresponds to variable J introduced in equation (14) in the report
     
    bigj = zeros(M,M*p); 
    bigj(1:M,1:M) = eye(M);
     
     
    %% ================ START GIBBS SAMPLING ==================================
     
    tic; % This is just a timer
    disp('Number of iterations');
     
     
    for irep = 1:nrep + nburn    % 7000 gibbs iterations starts here
        % Print iterations - this just updates on the progress of the sampling
        if mod(irep,it_print) == 0
            disp(irep);toc;
        end
        
    %% Draw 1: B_t from p(B_t|y,Sigma)
        
        % We use the function 'carter_kohn_hom' to to run the FFBS algorithm. 
        % This results in a 21x173 matrix, corresponding to one Gibbs sample
        % draw of each of the coefficients in each time period. The inputs
        % Sigmadraw and Qdraw are updated for each Gibbs sample repetition.
        
        [Btdraw] = carter_kohn_hom(y,Z,Sigmadraw,Qdraw,K,M,t,B_0_prmean,B_0_prvar);
        
      
    %% Draw 2: Q from p(Q^{-1}|y,B_t) which is i-Wishart
     
        % We draw Q from an Inverse Wishart distribution. The parameters 
        % of the distribution are derived as equation (11) in the main report.
        % The mean is taken as the inverse of the accumulated sum of squared 
        % errors added to the prior mean, and the variance is simply t.  
        
        % Differencing Btdraw to create the sum of squared errors
        
        Btemp = Btdraw(:,2:t)' - Btdraw(:,1:t-1)'; 
        sse_2Q = zeros(K,K);
        for i = 1:t-1
            sse_2Q = sse_2Q + Btemp(i,:)'*Btemp(i,:);
        end
     
        Qinv = inv(sse_2Q + Q_prmean);      % compute mean to use for Wishart draw
        Qinvdraw = wish(Qinv,t+Q_prvar);    % draw inv q from the wishart distribution
        Qdraw = inv(Qinvdraw);              % find non-inverse q 
        
    %% Draw 3: Sigma from p(Sigma|y,B_t) which is i-Wishart
     
        % We draw Sigma from an Inverse Wishart distribution. The parameters 
        % of the distirbution are derived as equation (10) in the main report. 
        % The mean is taken as the inverse of the sum of squared residuals
        % added to the prior mean. The variance is simply t. 
        
        % Find residuals using data and the current draw of coefficients
        
        resids = zeros(M,t);
        for i = 1:t
            resids(:,i) = y(:,i) - Z((i-1)*M+1:i*M,:)*Btdraw(:,i);
        end
        
        % Create a matrix for the accumulated sum of squared residuals, to
        % be used as the mean parameter in the i-Wishart draw below. 
        sse_2S = zeros(M,M);
        for i = 1:t
            sse_2S = sse_2S + resids(:,i)*resids(:,i)';
        end
        
        Sigmainv = inv(sse_2S + Sigma_prmean);          % compute mean to use for the Wishart
        Sigmainvdraw = wish(Sigmainv,t+Sigma_prvar);    % draw from the Wishsart distribution
        Sigmadraw = inv(Sigmainvdraw);                  % turn into non-inverse Sigma
        
    %% IRF 
        % We only apply IRF analysis once we have exceeded the burn-in draws phase.
        if irep > nburn;         
                
                % Create matrix that is going to contain all beta draws over
                % which we will take the mean after the Gibbs sampler as our moment
                % estimate: 
                Bt_postmean = Bt_postmean + Btdraw;
                
                % biga is the A matrix of the VAR(1) version of our VAR(2) model, 
                % found in equation (12. biga changes in every period of the 
                % analysis, because the coefficients are time varying, so we 
                % apply the analysis below in every time period. 
                
                biga = zeros(M*p,M*p); 
                for j = 1:p-1
                    biga(j*M+1:M*(j+1),M*(j-1)+1:j*M) = eye(M); % fill the A matrix with identity matrix (3) in bottom left corner
                end
     
                % The following procedure is applied separately in each time period. 
                
                % This loop takes coefficients of the relevant time period from
                % Bt_draw (which contains all coefficients for all t) and uses
                % them to update the biga matrix, so that it can change for
                % every t. 
                
                for i = 1:t 
                    bbtemp = Btdraw(M+1:K,i);  % get the draw of B(t) at time i=1,...,T  (exclude intercept)
                    splace = 0;
                    for ii = 1:p
                        for iii = 1:M
                            biga(iii,(ii-1)*M+1:ii*M) = bbtemp(splace+1:splace+M,1)'; % Load non-intercept coefficient draws
                            splace = splace + M;
                        end
                    end
                    
                    % Next we want to create a shock matrix in which the third
                    % column is [0 0 1]', therefore implementing a unit shock
                    % in the interest rate. 
                    
                    shock = eye(3);
                    
                    % Now get impulse responses for 1 through nhor future
                    % periods. impresp is a 3x63 matrix which contains 9
                    % response values in total for each period, 3 for each 
                    % variable. These three responses correspond to the 3
                    % possible shocks that are contained in the schock
                    % matrix
                    
                    % bigai is updated through mulitiplication with the 
                    % coefficient matrix after each time period. 
                                    
                    % Create a results matrix to store impulse responses in all periods
                    impresp = zeros(M,M*nhor); 
                    
                    % Fill in the first period of the results matrix with the shock (as defined above) 
                    impresp(1:M,1:M) = shock;
                    
                    % Create a separate variable for the a matrix so that we
                    % can update it for each period of the IRF analysis. 
                    bigai = biga; 
                    
                    % This follows the impulse response function as in equation 15.
                    % Fill in each period of the results matrix according to
                    % the impulse response function formula. 
                    
                    for j = 1:nhor-1
                        impresp(:,j*M+1:(j+1)*M) = bigj*bigai*bigj'*shock;
                        bigai = bigai*biga; % update the coefficient matrix for next period
                    end
     
                    % The section below keeps only the responses that we are interested in:
                    % - those from the periods 1975q1, 1981q3, and 1996q1
                    % - those that correspond to the shock in the interest 
                    % rate (i.e. those caused by the third column of our shock
                    % matrix). 
                    
                    if yearlab(i,1) == 1975.00;   % store only IRF from 1975:Q1
                        impf_m = zeros(M,nhor);
                        jj=0;
                        for ij = 1:nhor
                            jj = jj + M;    % select only the third column for each time period of the IRF
                            impf_m(:,ij) = impresp(:,jj);
                        end           
                    % For each iteration of the Gibbs sampler, fill in the
                    % results along the first dimension 
                        imp75(irep-nburn,:,:) = impf_m; 
                    end
                    
                    if yearlab(i,1) == 1981.50;   % store only IRF from 1981:Q3
                        impf_m = zeros(M,nhor);
                        jj=0;
                        for ij = 1:nhor
                            jj = jj + M;    % select only the third column for each time period of the IRF
                            impf_m(:,ij) = impresp(:,jj);
                        end
                    % For each iteration of the Gibbs sample, fill in the
                    % results along the first dimension 
                        imp81(irep-nburn,:,:) = impf_m; 
                    end
                    
                    if yearlab(i,1) == 1996.00;   % store only IRF from 1996:Q1
                        impf_m = zeros(M,nhor);
                        jj=0;
                        for ij = 1:nhor
                            jj = jj + M;    % select only the third column for each time period of the IRF
                            impf_m(:,ij) = impresp(:,jj);
                        end
                    % For each iteration of the Gibbs sample, fill in the
                    % results along the first dimension 
                        imp96(irep-nburn,:,:) = impf_m; 
                    end
                end % End getting impulses for each time period 
            end % End the impulse response calculation section   
    end % End main Gibbs loop (for irep = 1:nrep+nburn)
     
    clc;
    toc; % Stop timer and print total time
    %% ================ END GIBBS SAMPLING ==================================
     
    % Even though it is not used in our IRF analysis since we are integrating
    % that into the Gibbs sampling loop, here is how to take the mean of the 
    % draw of the betas as moment estimate:
     
    Bt_postmean = Bt_postmean./nrep;
     
    %% Graphs and tables
     
    % This works out the percentage range of that each variable's coefficient
    % spans across time. I.e. to what extent was the TVP facility used by each
    % variable in the model? This is calculated by finding the range for each
    % variable as a percentage of the mean coefficient size for that variable.
    % The result is a 21x1 vector, and it is reshaped into a 3x7 matrix in order
    % to map onto the system of equations (2,3,and 4) set out in the report.
     
    Bt_range = ones(21,1)
    for i = 1:21
        Bt_range(i) = abs((max(Bt_postmean(i,:))-min(Bt_postmean(i,:)))/mean(Bt_postmean(i,:)))
    end 
    Bt_range = reshape(Bt_range,3,7)
     
    % Create a table of coefficient ranges for export to the report
     
    rowNames = {'Inflation','Unemployment','Interest Rate'};
    colNames = {'Intercept','Inf_1','Unemp_1', 'IR_1','Inf_2','Unemp_2', 'IR_2'};
    pc_change_table = array2table(Bt_range,'RowNames',rowNames,'VariableNames',colNames)
     
    writetable(pc_change_table,'pc_change.csv')
     
    % Now plot a separate chart for each of the coefficients
     
    figure
    for i = 1:21       
        subplot(7,3,i)
        plot(1:t,Bt_postmean(i,:))  
    end 
     
    % Now we move to plotting the IRF. This section takes moments along the
    % first dimension, i.e. across the Gibbs sample iterations. The moments
    % are for the 16th, 50th and 84th percentile. 
     
        qus = [.16, .5, .84];
        imp75XY=squeeze(quantile(imp75,qus)); 
        imp81XY=squeeze(quantile(imp81,qus));
        imp96XY=squeeze(quantile(imp96,qus));
        
        % Plot impulse responses
        figure       
        set(0,'DefaultAxesColorOrder',[0 0 0],...
            'DefaultAxesLineStyleOrder','--|-|--')
        subplot(3,3,1)
        plot(1:nhor,squeeze(imp75XY(:,1,:)))
        title('Impulse response of inflation, 1975:Q1')
        xlim([1 nhor])
        ylim([-0.2 0.1])
    %    % yline(0)
        set(gca,'XTick',0:3:nhor)
        subplot(3,3,2)
        plot(1:nhor,squeeze(imp75XY(:,2,:)))
        title('Impulse response of unemployment, 1975:Q1')
        xlim([1 nhor])
        ylim([-0.2 0.2])
       % yline(0)
        set(gca,'XTick',0:3:nhor)    
        subplot(3,3,3)
        %ylim([0 1])
    %    % yline(0)
        plot(1:nhor,squeeze(imp75XY(:,3,:)))
        title('Impulse response of interest rate, 1975:Q1')
        xlim([1 nhor])
        %ylim([-0.3 0.1])
       % yline(0)
        set(gca,'XTick',0:3:nhor)    
        subplot(3,3,4)
        plot(1:nhor,squeeze(imp81XY(:,1,:)))
        title('Impulse response of inflation, 1981:Q3')
        xlim([1 nhor])
        ylim([-0.2 0.1])
       % yline(0)
        set(gca,'XTick',0:3:nhor)    
        subplot(3,3,5)
        plot(1:nhor,squeeze(imp81XY(:,2,:)))
        title('Impulse response of unemployment, 1981:Q3')
        xlim([1 nhor])
        ylim([-0.2 0.2])
       % yline(0)
        set(gca,'XTick',0:3:nhor)    
        subplot(3,3,6)
        plot(1:nhor,squeeze(imp81XY(:,3,:)))
        title('Impulse response of interest rate, 1981:Q3')
        xlim([1 nhor])
        %ylim([-0.4 0.1])
       % yline(0)
        set(gca,'XTick',0:3:nhor)    
        subplot(3,3,7)
        plot(1:nhor,squeeze(imp96XY(:,1,:)))
        title('Impulse response of inflation, 1996:Q1')
        xlim([1 nhor])
        ylim([-0.2 0.1])
       % yline(0)
        set(gca,'XTick',0:3:nhor)    
        subplot(3,3,8)
        plot(1:nhor,squeeze(imp96XY(:,2,:)))
        title('Impulse response of unemployment, 1996:Q1')
        xlim([1 nhor])
        ylim([-0.2 0.2])
       % yline(0)
        set(gca,'XTick',0:3:nhor)
        subplot(3,3,9)
        plot(1:nhor,squeeze(imp96XY(:,3,:)))
        title('Impulse response of interest rate, 1996:Q1')
        xlim([1 nhor])
         %ylim([0 1])
        % yline(0)
        set(gca,'XTick',0:3:nhor)
        
    disp('             ')
    disp('To plot impulse responses, use:         plot(1:nhor,squeeze(imp75XY(:,VAR,:)))           ')
    disp('             ')
    disp('where VAR=1 for impulses of inflation, VAR=2 for unemployment and VAR=3 for interest rate')
    

    三、运行结果

    在这里插入图片描述
    在这里插入图片描述

    四、备注

    版本:2014a

    展开全文
  • VSC-HVDC的matlab模型

    2018-10-30 09:11:38
    描述 200 MVA(+/- 100 kV DC)强制换向电压源转换器(VSC)互连用于将功率从230 kV,2000 MVA,50 Hz系统传输到另一个相同的AC...请注意,模型的“模型初始化”功能会自动在MATLAB®工作空间中设置这两个采样时间。
  • matlab开发-反转var参数Intomap参数。将向量自回归模型参数转化为移动平均模型参数的函数。
  • nmf的matlab代码 ...由于var模型的数据是时间序列的,因此数据分布需要创建块,因此可以对问题进行块随机化。 在var_vectorize.cpp文件中创建一个分布式矢量化模块。 分布式kronecker产品位于var_kr
  • 一、简介 VAR:向量自回归模型,结果仅具有统计上的 意义 SVAR:结构向量自回归模型 TVP-VAR:Time Varying Parameter-Stochastic Volatility-Vector Auto Regression。...在VAR模型中,所有的参数遵循一阶游走
  • 一、简介 VAR:向量自回归模型,结果仅具有统计上的 意义 SVAR:结构向量自回归模型 TVP-VAR:Time Varying Parameter-Stochastic Volatility-Vector Auto Regression。...在VAR模型中,所有的参数遵循一阶游走
  • (S)VAR/协整模型的类型直观明显时,这通常很有用。 第 3 节对单位根进行 ADF 测试(正在进行中)。 第 4 节对常规 VAR 进行测试,以选择最佳腿长。 第 5 部分使用经济直觉和判断力选择顶部的滞后长度,然后为您绘制/...
  • ar模型matlab代码MATLAB中的11种经典时间序列预测方法 凯文 机器学习实施的兴起,引起了不同行业的兴趣,将其用于时间序列问题的分类和预测。...VAR模型预测马来西亚/美国的汇率 使用ARIMA进行库存预测 使用ARIMA和
  • 致谢:我们依靠惩罚的LP代码,以及依靠VAR模型求平均的代码。 我们对这两组代码进行了一些修改,以改善其运行时间,而不会影响其数字输出。 我们还使用的动态因子模型代码和数据。 内容 :纸,补充资料和文档 :主纸...
  • MATLAB用拟合出的代码绘图自述文件 这组GAUSS代码可计算由真实商业周期(RBC)模型生成...DGP(在存档中保存为Simulation.txt)的200条观测值的样本路径,然后拟合有限滞后VAR模型并计算对结构冲击的冲激响应使用Blanc
  • matlab自相关代码真棒VAR 向量自回归资源的精选列表。 目录: 的MATLAB 工具箱 :向量自回归模型 :收集Matlab例程以执行VAR分析(Ambrogio Cesa-Bianchi) :经验宏工具箱(F. Ferroni和F. Canova) :宏观经济建模...

空空如也

空空如也

1 2 3 4 5 6
收藏数 117
精华内容 46
关键字:

matlabvar模型

matlab 订阅