%==================================4.1 Main================================
clear;
clc;
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
Matlab TVP-VAR-project时变参数随机波动率向量自回归模型.zip
版权申诉
5星 · 超过95%的资源 48 浏览量
2021-12-26
12:55:28
上传
评论 4
收藏 38KB ZIP 举报
天天Matlab科研工作室
- 粉丝: 2w+
- 资源: 7248
最新资源
- 电机控制霍尔传感器和反电动势的关系分析
- 计算机生产实习:OA后台管理(web前端+Java后端)压缩文件包
- greenplum6.10安装时缺失的fpm包
- 缺陷检测-轻量化PCB表面缺陷检测算法实现-工业视觉+六大缺陷检出-优质项目实战.zip
- JAVA 中的Spring框架介绍包括起源、体系结构、核心部分、特点等
- 2024年小米汽车产业链分析及新品上市全景洞察报告
- 基于Qt和C++实现的偏3D风格的异型窗体界面操作+源码(期末大作业&课设&项目开发)
- 基于yolov8的花卉分类系统,包含训练好的权重和推理代码,GUI界面,支持图片、视频、摄像头输入,支持检测结果导出
- 基于图形化编程的单片机教学案例研究
- 基于matlab语音识别的信号灯图像模拟控制技术代码19
资源上传下载、课程学习等过程中有任何疑问或建议,欢迎提出宝贵意见哦~我们会及时处理!
点击此处反馈
- 1
- 2
- 3
前往页