R 中的历史分解
Historical Decomposition In R
我目前正在尝试 运行 在 R 中对我的数据系列进行历史分解。
我阅读了大量论文，它们都对如何进行历史分解提供了以下解释：
其中右侧的总和是 Yt+k 的 "dynamic forecast" 或 "base projection"，条件是在时间 t 可用的信息。左边的总和是实际序列与基础预测之间的差异，由于 t+1 到 t+k 期间变量的创新
我对基本投影感到非常困惑，我不确定使用的是什么数据！
我的尝试。
我有一个 6 变量 VAR，有 55 个观察值。我使用 Cholesky 分解得到模型的结构形式。完成此操作后，我使用 Phi 函数来获取 SVAR 的结构移动平均表示。然后我存储这个 Phi "array" 以便我以后可以使用它。
varFT <- VAR(Enddata[,c(2,3,4,5,6,7)], p = 4, type = c("const"))
Amat <- diag(6)
Amat
Bmat <- diag(6)
Bmat[1,1] <- NA
Bmat[2,2] <- NA
Bmat[3,3] <- NA
Bmat[4,4] <- NA
Bmat[5,5] <- NA
Bmat[6,6] <- NA
#play around with col/row names to make them pretty/understandable.
colnames(Bmat) <- c("G", "FT", "T","R", "P", "Y")
rownames(Bmat) <- c("G", "FT", "T", "R", "P", "Y")
Amat[1,5] <- 0
Amat[1,4] <- 0
Amat[1,3] <- 0
#Make Amat lower triangular, leave Bmat as diag.
Amat[5,1:4] <- NA
Amat[4, 1:3] <- NA
Amat[3,1:2] <- NA
Amat[2,1] <- NA
Amat[6,1:5] <- NA
svarFT <- SVAR(varFT, estmethod = c("scoring"), Amat = Amat, Bmat = Bmat)
MA <- Phi(svarFT, nstep = 55)
MAarray <- function(x){
resid_store = array(0, dim=c(6,6,54))
resid_store[,,1] = (Phi(x, nstep = 54))[,,1]
for (d in 1:54){
resid_store[,,d] = Phi(x,nstep = 54)[,,d]
}
return(resid_store)
}
Part1 <-MAarray(MA)
我认为我已经获得了基本投影所需的信息，但我不知道从这里到哪里去。
目标
我想做的是评估 VAR 中的第一个变量在整个样本期间对 VAR 中的第六个变量的影响。
如有任何帮助，我们将不胜感激。
历史分解真正解决了一个系列的 错误 如何影响 VAR 中的另一个系列。最简单的方法是创建一个拟合误差数组。从这里开始，您将需要一个三层嵌套的 for 循环：
循环拟合冲击系列：for (iShock in 1:6)
循环给定拟合冲击的时间维度，从基期之后的周期开始：for (iShockPeriod in 1:55)
模拟个体实现该冲击值对其余样本的影响：for (iResponsePeriod in iShockPeriod:55)
您实际上得到了一个尺寸为（例如）6x6x55x55 的 4 维数组。 (i,j,k,l) 元素类似于第 k 期第 i 个系列的冲击对第 l 期第 j 个系列的影响。当我以前写过实现时，在你进行的过程中至少对这些维度进行求和通常是有意义的，这样你就不会得到这么大的数组。
不幸的是，我没有 R 中的实现可以分享，但我一直在 Stata 中开发一个。如果我很快就可以使用它，我会用 link 更新它。
我将函数 VARhd 从 Cesa-Bianchi's Matlab Toolbox 翻译成 R 代码。我的函数与 R.
中 vars 包中的 VAR 函数兼容
MATLAB中的原函数：
function HD = VARhd(VAR,VARopt)
% =======================================================================
% Computes the historical decomposition of the times series in a VAR
% estimated with VARmodel and identified with VARir/VARfevd
% =======================================================================
% HD = VARhd(VAR,VARopt)
% -----------------------------------------------------------------------
% INPUTS
% - VAR: VAR results obtained with VARmodel (structure)
% - VARopt: options of the IRFs (see VARoption)
% OUTPUT
% - HD(t,j,k): matrix with 't' steps, containing the IRF of 'j' variable
% to 'k' shock
% - VARopt: options of the IRFs (see VARoption)
% =======================================================================
% Ambrogio Cesa Bianchi, April 2014
% ambrogio.cesabianchi@gmail.com
%% Check inputs
%===============================================
if ~exist('VARopt','var')
error('You need to provide VAR options (VARopt from VARmodel)');
end
%% Retrieve and initialize variables
%=============================================================
invA = VARopt.invA; % inverse of the A matrix
Fcomp = VARopt.Fcomp; % Companion matrix
det = VAR.det; % constant and/or trends
F = VAR.Ft'; % make comparable to notes
eps = invA\transpose(VAR.residuals); % structural errors
nvar = VAR.nvar; % number of endogenous variables
nvarXeq = VAR.nvar * VAR.nlag; % number of lagged endogenous per equation
nlag = VAR.nlag; % number of lags
nvar_ex = VAR.nvar_ex; % number of exogenous (excluding constant and trend)
Y = VAR.Y; % left-hand side
X = VAR.X(:,1+det:nvarXeq+det); % right-hand side (no exogenous)
nobs = size(Y,1); % number of observations
%% Compute historical decompositions
%===================================
% Contribution of each shock
invA_big = zeros(nvarXeq,nvar);
invA_big(1:nvar,:) = invA;
Icomp = [eye(nvar) zeros(nvar,(nlag-1)*nvar)];
HDshock_big = zeros(nlag*nvar,nobs+1,nvar);
HDshock = zeros(nvar,nobs+1,nvar);
for j=1:nvar; % for each variable
eps_big = zeros(nvar,nobs+1); % matrix of shocks conformable with companion
eps_big(j,2:end) = eps(j,:);
for i = 2:nobs+1
HDshock_big(:,i,j) = invA_big*eps_big(:,i) + Fcomp*HDshock_big(:,i-1,j);
HDshock(:,i,j) = Icomp*HDshock_big(:,i,j);
end
end
% Initial value
HDinit_big = zeros(nlag*nvar,nobs+1);
HDinit = zeros(nvar, nobs+1);
HDinit_big(:,1) = X(1,:)';
HDinit(:,1) = Icomp*HDinit_big(:,1);
for i = 2:nobs+1
HDinit_big(:,i) = Fcomp*HDinit_big(:,i-1);
HDinit(:,i) = Icomp *HDinit_big(:,i);
end
% Constant
HDconst_big = zeros(nlag*nvar,nobs+1);
HDconst = zeros(nvar, nobs+1);
CC = zeros(nlag*nvar,1);
if det>0
CC(1:nvar,:) = F(:,1);
for i = 2:nobs+1
HDconst_big(:,i) = CC + Fcomp*HDconst_big(:,i-1);
HDconst(:,i) = Icomp * HDconst_big(:,i);
end
end
% Linear trend
HDtrend_big = zeros(nlag*nvar,nobs+1);
HDtrend = zeros(nvar, nobs+1);
TT = zeros(nlag*nvar,1);
if det>1;
TT(1:nvar,:) = F(:,2);
for i = 2:nobs+1
HDtrend_big(:,i) = TT*(i-1) + Fcomp*HDtrend_big(:,i-1);
HDtrend(:,i) = Icomp * HDtrend_big(:,i);
end
end
% Quadratic trend
HDtrend2_big = zeros(nlag*nvar, nobs+1);
HDtrend2 = zeros(nvar, nobs+1);
TT2 = zeros(nlag*nvar,1);
if det>2;
TT2(1:nvar,:) = F(:,3);
for i = 2:nobs+1
HDtrend2_big(:,i) = TT2*((i-1)^2) + Fcomp*HDtrend2_big(:,i-1);
HDtrend2(:,i) = Icomp * HDtrend2_big(:,i);
end
end
% Exogenous
HDexo_big = zeros(nlag*nvar,nobs+1);
HDexo = zeros(nvar,nobs+1);
EXO = zeros(nlag*nvar,nvar_ex);
if nvar_ex>0;
VARexo = VAR.X_EX;
EXO(1:nvar,:) = F(:,nvar*nlag+det+1:end); % this is c in my notes
for i = 2:nobs+1
HDexo_big(:,i) = EXO*VARexo(i-1,:)' + Fcomp*HDexo_big(:,i-1);
HDexo(:,i) = Icomp * HDexo_big(:,i);
end
end
% All decompositions must add up to the original data
HDendo = HDinit + HDconst + HDtrend + HDtrend2 + HDexo + sum(HDshock,3);
%% Save and reshape all HDs
%==========================
HD.shock = zeros(nobs+nlag,nvar,nvar); % [nobs x shock x var]
for i=1:nvar
for j=1:nvar
HD.shock(:,j,i) = [nan(nlag,1); HDshock(i,2:end,j)'];
end
end
HD.init = [nan(nlag-1,nvar); HDinit(:,1:end)']; % [nobs x var]
HD.const = [nan(nlag,nvar); HDconst(:,2:end)']; % [nobs x var]
HD.trend = [nan(nlag,nvar); HDtrend(:,2:end)']; % [nobs x var]
HD.trend2 = [nan(nlag,nvar); HDtrend2(:,2:end)']; % [nobs x var]
HD.exo = [nan(nlag,nvar); HDexo(:,2:end)']; % [nobs x var]
HD.endo = [nan(nlag,nvar); HDendo(:,2:end)']; % [nobs x var]
我的 R 版本（基于 vars 包）：
VARhd <- function(Estimation){
## make X and Y
nlag <- Estimation$p # number of lags
DATA <- Estimation$y # data
QQ <- VARmakexy(DATA,nlag,1)
## Retrieve and initialize variables
invA <- t(chol(as.matrix(summary(Estimation)$covres))) # inverse of the A matrix
Fcomp <- companionmatrix(Estimation) # Companion matrix
#det <- c_case # constant and/or trends
F1 <- t(QQ$Ft) # make comparable to notes
eps <- ginv(invA) %*% t(residuals(Estimation)) # structural errors
nvar <- Estimation$K # number of endogenous variables
nvarXeq <- nvar * nlag # number of lagged endogenous per equation
nvar_ex <- 0 # number of exogenous (excluding constant and trend)
Y <- QQ$Y # left-hand side
#X <- QQ$X[,(1+det):(nvarXeq+det)] # right-hand side (no exogenous)
nobs <- nrow(Y) # number of observations
## Compute historical decompositions
# Contribution of each shock
invA_big <- matrix(0,nvarXeq,nvar)
invA_big[1:nvar,] <- invA
Icomp <- cbind(diag(nvar), matrix(0,nvar,(nlag-1)*nvar))
HDshock_big <- array(0, dim=c(nlag*nvar,nobs+1,nvar))
HDshock <- array(0, dim=c(nvar,(nobs+1),nvar))
for (j in 1:nvar){ # for each variable
eps_big <- matrix(0,nvar,(nobs+1)) # matrix of shocks conformable with companion
eps_big[j,2:ncol(eps_big)] <- eps[j,]
for (i in 2:(nobs+1)){
HDshock_big[,i,j] <- invA_big %*% eps_big[,i] + Fcomp %*% HDshock_big[,(i-1),j]
HDshock[,i,j] <- Icomp %*% HDshock_big[,i,j]
}
}
HD.shock <- array(0, dim=c((nobs+nlag),nvar,nvar)) # [nobs x shock x var]
for (i in 1:nvar){
for (j in 1:nvar){
HD.shock[,j,i] <- c(rep(NA,nlag), HDshock[i,(2:dim(HDshock)[2]),j])
}
}
return(HD.shock)
}
作为输入参数，您必须使用 R 中 vars 包中的 VAR 函数。函数 returns 一个 3 维数组：观察数 x 冲击数 x 变量数。
（注意：我没有翻译整个函数，例如我省略了外生变量的情况。）
要 运行 它，您需要两个额外的函数，这两个函数也是从 Bianchi 的工具箱翻译而来的：
VARmakexy <- function(DATA,lags,c_case){
nobs <- nrow(DATA)
#Y matrix
Y <- DATA[(lags+1):nrow(DATA),]
Y <- DATA[-c(1:lags),]
#X-matrix
if (c_case==0){
X <- NA
for (jj in 0:(lags-1)){
X <- rbind(DATA[(jj+1):(nobs-lags+jj),])
}
} else if(c_case==1){ #constant
X <- NA
for (jj in 0:(lags-1)){
X <- rbind(DATA[(jj+1):(nobs-lags+jj),])
}
X <- cbind(matrix(1,(nobs-lags),1), X)
} else if(c_case==2){ # time trend and constant
X <- NA
for (jj in 0:(lags-1)){
X <- rbind(DATA[(jj+1):(nobs-lags+jj),])
}
trend <- c(1:nrow(X))
X <-cbind(matrix(1,(nobs-lags),1), t(trend))
}
A <- (t(X) %*% as.matrix(X))
B <- (as.matrix(t(X)) %*% as.matrix(Y))
Ft <- ginv(A) %*% B
retu <- list(X=X,Y=Y, Ft=Ft)
return(retu)
}
companionmatrix <- function (x)
{
if (!(class(x) == "varest")) {
stop("\nPlease provide an object of class 'varest', generated by 'VAR()'.\n")
}
K <- x$K
p <- x$p
A <- unlist(Acoef(x))
companion <- matrix(0, nrow = K * p, ncol = K * p)
companion[1:K, 1:(K * p)] <- A
if (p > 1) {
j <- 0
for (i in (K + 1):(K * p)) {
j <- j + 1
companion[i, j] <- 1
}
}
return(companion)
}
这是一个简短的例子：
library(vars)
data(Canada)
ab<-VAR(Canada, p = 2, type = "both")
HD <- VARhd(Estimation=ab)
HD[,,1] # historical decomposition of the first variable (employment)
这是excel中的情节：
我目前正在尝试 运行 在 R 中对我的数据系列进行历史分解。
我阅读了大量论文，它们都对如何进行历史分解提供了以下解释：
其中右侧的总和是 Yt+k 的 "dynamic forecast" 或 "base projection"，条件是在时间 t 可用的信息。左边的总和是实际序列与基础预测之间的差异，由于 t+1 到 t+k 期间变量的创新
我对基本投影感到非常困惑，我不确定使用的是什么数据！
我的尝试。
我有一个 6 变量 VAR，有 55 个观察值。我使用 Cholesky 分解得到模型的结构形式。完成此操作后，我使用 Phi 函数来获取 SVAR 的结构移动平均表示。然后我存储这个 Phi "array" 以便我以后可以使用它。
varFT <- VAR(Enddata[,c(2,3,4,5,6,7)], p = 4, type = c("const"))
Amat <- diag(6)
Amat
Bmat <- diag(6)
Bmat[1,1] <- NA
Bmat[2,2] <- NA
Bmat[3,3] <- NA
Bmat[4,4] <- NA
Bmat[5,5] <- NA
Bmat[6,6] <- NA
#play around with col/row names to make them pretty/understandable.
colnames(Bmat) <- c("G", "FT", "T","R", "P", "Y")
rownames(Bmat) <- c("G", "FT", "T", "R", "P", "Y")
Amat[1,5] <- 0
Amat[1,4] <- 0
Amat[1,3] <- 0
#Make Amat lower triangular, leave Bmat as diag.
Amat[5,1:4] <- NA
Amat[4, 1:3] <- NA
Amat[3,1:2] <- NA
Amat[2,1] <- NA
Amat[6,1:5] <- NA
svarFT <- SVAR(varFT, estmethod = c("scoring"), Amat = Amat, Bmat = Bmat)
MA <- Phi(svarFT, nstep = 55)
MAarray <- function(x){
resid_store = array(0, dim=c(6,6,54))
resid_store[,,1] = (Phi(x, nstep = 54))[,,1]
for (d in 1:54){
resid_store[,,d] = Phi(x,nstep = 54)[,,d]
}
return(resid_store)
}
Part1 <-MAarray(MA)
我认为我已经获得了基本投影所需的信息，但我不知道从这里到哪里去。
目标 我想做的是评估 VAR 中的第一个变量在整个样本期间对 VAR 中的第六个变量的影响。
如有任何帮助，我们将不胜感激。
历史分解真正解决了一个系列的 错误 如何影响 VAR 中的另一个系列。最简单的方法是创建一个拟合误差数组。从这里开始，您将需要一个三层嵌套的 for 循环：
循环拟合冲击系列：
for (iShock in 1:6)循环给定拟合冲击的时间维度，从基期之后的周期开始：
for (iShockPeriod in 1:55)模拟个体实现该冲击值对其余样本的影响：
for (iResponsePeriod in iShockPeriod:55)
您实际上得到了一个尺寸为（例如）6x6x55x55 的 4 维数组。 (i,j,k,l) 元素类似于第 k 期第 i 个系列的冲击对第 l 期第 j 个系列的影响。当我以前写过实现时，在你进行的过程中至少对这些维度进行求和通常是有意义的，这样你就不会得到这么大的数组。
不幸的是，我没有 R 中的实现可以分享，但我一直在 Stata 中开发一个。如果我很快就可以使用它，我会用 link 更新它。
我将函数 VARhd 从 Cesa-Bianchi's Matlab Toolbox 翻译成 R 代码。我的函数与 R.
vars 包中的 VAR 函数兼容
MATLAB中的原函数：
function HD = VARhd(VAR,VARopt)
% =======================================================================
% Computes the historical decomposition of the times series in a VAR
% estimated with VARmodel and identified with VARir/VARfevd
% =======================================================================
% HD = VARhd(VAR,VARopt)
% -----------------------------------------------------------------------
% INPUTS
% - VAR: VAR results obtained with VARmodel (structure)
% - VARopt: options of the IRFs (see VARoption)
% OUTPUT
% - HD(t,j,k): matrix with 't' steps, containing the IRF of 'j' variable
% to 'k' shock
% - VARopt: options of the IRFs (see VARoption)
% =======================================================================
% Ambrogio Cesa Bianchi, April 2014
% ambrogio.cesabianchi@gmail.com
%% Check inputs
%===============================================
if ~exist('VARopt','var')
error('You need to provide VAR options (VARopt from VARmodel)');
end
%% Retrieve and initialize variables
%=============================================================
invA = VARopt.invA; % inverse of the A matrix
Fcomp = VARopt.Fcomp; % Companion matrix
det = VAR.det; % constant and/or trends
F = VAR.Ft'; % make comparable to notes
eps = invA\transpose(VAR.residuals); % structural errors
nvar = VAR.nvar; % number of endogenous variables
nvarXeq = VAR.nvar * VAR.nlag; % number of lagged endogenous per equation
nlag = VAR.nlag; % number of lags
nvar_ex = VAR.nvar_ex; % number of exogenous (excluding constant and trend)
Y = VAR.Y; % left-hand side
X = VAR.X(:,1+det:nvarXeq+det); % right-hand side (no exogenous)
nobs = size(Y,1); % number of observations
%% Compute historical decompositions
%===================================
% Contribution of each shock
invA_big = zeros(nvarXeq,nvar);
invA_big(1:nvar,:) = invA;
Icomp = [eye(nvar) zeros(nvar,(nlag-1)*nvar)];
HDshock_big = zeros(nlag*nvar,nobs+1,nvar);
HDshock = zeros(nvar,nobs+1,nvar);
for j=1:nvar; % for each variable
eps_big = zeros(nvar,nobs+1); % matrix of shocks conformable with companion
eps_big(j,2:end) = eps(j,:);
for i = 2:nobs+1
HDshock_big(:,i,j) = invA_big*eps_big(:,i) + Fcomp*HDshock_big(:,i-1,j);
HDshock(:,i,j) = Icomp*HDshock_big(:,i,j);
end
end
% Initial value
HDinit_big = zeros(nlag*nvar,nobs+1);
HDinit = zeros(nvar, nobs+1);
HDinit_big(:,1) = X(1,:)';
HDinit(:,1) = Icomp*HDinit_big(:,1);
for i = 2:nobs+1
HDinit_big(:,i) = Fcomp*HDinit_big(:,i-1);
HDinit(:,i) = Icomp *HDinit_big(:,i);
end
% Constant
HDconst_big = zeros(nlag*nvar,nobs+1);
HDconst = zeros(nvar, nobs+1);
CC = zeros(nlag*nvar,1);
if det>0
CC(1:nvar,:) = F(:,1);
for i = 2:nobs+1
HDconst_big(:,i) = CC + Fcomp*HDconst_big(:,i-1);
HDconst(:,i) = Icomp * HDconst_big(:,i);
end
end
% Linear trend
HDtrend_big = zeros(nlag*nvar,nobs+1);
HDtrend = zeros(nvar, nobs+1);
TT = zeros(nlag*nvar,1);
if det>1;
TT(1:nvar,:) = F(:,2);
for i = 2:nobs+1
HDtrend_big(:,i) = TT*(i-1) + Fcomp*HDtrend_big(:,i-1);
HDtrend(:,i) = Icomp * HDtrend_big(:,i);
end
end
% Quadratic trend
HDtrend2_big = zeros(nlag*nvar, nobs+1);
HDtrend2 = zeros(nvar, nobs+1);
TT2 = zeros(nlag*nvar,1);
if det>2;
TT2(1:nvar,:) = F(:,3);
for i = 2:nobs+1
HDtrend2_big(:,i) = TT2*((i-1)^2) + Fcomp*HDtrend2_big(:,i-1);
HDtrend2(:,i) = Icomp * HDtrend2_big(:,i);
end
end
% Exogenous
HDexo_big = zeros(nlag*nvar,nobs+1);
HDexo = zeros(nvar,nobs+1);
EXO = zeros(nlag*nvar,nvar_ex);
if nvar_ex>0;
VARexo = VAR.X_EX;
EXO(1:nvar,:) = F(:,nvar*nlag+det+1:end); % this is c in my notes
for i = 2:nobs+1
HDexo_big(:,i) = EXO*VARexo(i-1,:)' + Fcomp*HDexo_big(:,i-1);
HDexo(:,i) = Icomp * HDexo_big(:,i);
end
end
% All decompositions must add up to the original data
HDendo = HDinit + HDconst + HDtrend + HDtrend2 + HDexo + sum(HDshock,3);
%% Save and reshape all HDs
%==========================
HD.shock = zeros(nobs+nlag,nvar,nvar); % [nobs x shock x var]
for i=1:nvar
for j=1:nvar
HD.shock(:,j,i) = [nan(nlag,1); HDshock(i,2:end,j)'];
end
end
HD.init = [nan(nlag-1,nvar); HDinit(:,1:end)']; % [nobs x var]
HD.const = [nan(nlag,nvar); HDconst(:,2:end)']; % [nobs x var]
HD.trend = [nan(nlag,nvar); HDtrend(:,2:end)']; % [nobs x var]
HD.trend2 = [nan(nlag,nvar); HDtrend2(:,2:end)']; % [nobs x var]
HD.exo = [nan(nlag,nvar); HDexo(:,2:end)']; % [nobs x var]
HD.endo = [nan(nlag,nvar); HDendo(:,2:end)']; % [nobs x var]
我的 R 版本（基于 vars 包）：
VARhd <- function(Estimation){
## make X and Y
nlag <- Estimation$p # number of lags
DATA <- Estimation$y # data
QQ <- VARmakexy(DATA,nlag,1)
## Retrieve and initialize variables
invA <- t(chol(as.matrix(summary(Estimation)$covres))) # inverse of the A matrix
Fcomp <- companionmatrix(Estimation) # Companion matrix
#det <- c_case # constant and/or trends
F1 <- t(QQ$Ft) # make comparable to notes
eps <- ginv(invA) %*% t(residuals(Estimation)) # structural errors
nvar <- Estimation$K # number of endogenous variables
nvarXeq <- nvar * nlag # number of lagged endogenous per equation
nvar_ex <- 0 # number of exogenous (excluding constant and trend)
Y <- QQ$Y # left-hand side
#X <- QQ$X[,(1+det):(nvarXeq+det)] # right-hand side (no exogenous)
nobs <- nrow(Y) # number of observations
## Compute historical decompositions
# Contribution of each shock
invA_big <- matrix(0,nvarXeq,nvar)
invA_big[1:nvar,] <- invA
Icomp <- cbind(diag(nvar), matrix(0,nvar,(nlag-1)*nvar))
HDshock_big <- array(0, dim=c(nlag*nvar,nobs+1,nvar))
HDshock <- array(0, dim=c(nvar,(nobs+1),nvar))
for (j in 1:nvar){ # for each variable
eps_big <- matrix(0,nvar,(nobs+1)) # matrix of shocks conformable with companion
eps_big[j,2:ncol(eps_big)] <- eps[j,]
for (i in 2:(nobs+1)){
HDshock_big[,i,j] <- invA_big %*% eps_big[,i] + Fcomp %*% HDshock_big[,(i-1),j]
HDshock[,i,j] <- Icomp %*% HDshock_big[,i,j]
}
}
HD.shock <- array(0, dim=c((nobs+nlag),nvar,nvar)) # [nobs x shock x var]
for (i in 1:nvar){
for (j in 1:nvar){
HD.shock[,j,i] <- c(rep(NA,nlag), HDshock[i,(2:dim(HDshock)[2]),j])
}
}
return(HD.shock)
}
作为输入参数，您必须使用 R 中 vars 包中的 VAR 函数。函数 returns 一个 3 维数组：观察数 x 冲击数 x 变量数。
（注意：我没有翻译整个函数，例如我省略了外生变量的情况。）
要 运行 它，您需要两个额外的函数，这两个函数也是从 Bianchi 的工具箱翻译而来的：
VARmakexy <- function(DATA,lags,c_case){
nobs <- nrow(DATA)
#Y matrix
Y <- DATA[(lags+1):nrow(DATA),]
Y <- DATA[-c(1:lags),]
#X-matrix
if (c_case==0){
X <- NA
for (jj in 0:(lags-1)){
X <- rbind(DATA[(jj+1):(nobs-lags+jj),])
}
} else if(c_case==1){ #constant
X <- NA
for (jj in 0:(lags-1)){
X <- rbind(DATA[(jj+1):(nobs-lags+jj),])
}
X <- cbind(matrix(1,(nobs-lags),1), X)
} else if(c_case==2){ # time trend and constant
X <- NA
for (jj in 0:(lags-1)){
X <- rbind(DATA[(jj+1):(nobs-lags+jj),])
}
trend <- c(1:nrow(X))
X <-cbind(matrix(1,(nobs-lags),1), t(trend))
}
A <- (t(X) %*% as.matrix(X))
B <- (as.matrix(t(X)) %*% as.matrix(Y))
Ft <- ginv(A) %*% B
retu <- list(X=X,Y=Y, Ft=Ft)
return(retu)
}
companionmatrix <- function (x)
{
if (!(class(x) == "varest")) {
stop("\nPlease provide an object of class 'varest', generated by 'VAR()'.\n")
}
K <- x$K
p <- x$p
A <- unlist(Acoef(x))
companion <- matrix(0, nrow = K * p, ncol = K * p)
companion[1:K, 1:(K * p)] <- A
if (p > 1) {
j <- 0
for (i in (K + 1):(K * p)) {
j <- j + 1
companion[i, j] <- 1
}
}
return(companion)
}
这是一个简短的例子：
library(vars)
data(Canada)
ab<-VAR(Canada, p = 2, type = "both")
HD <- VARhd(Estimation=ab)
HD[,,1] # historical decomposition of the first variable (employment)
这是excel中的情节：