Историческое разложение в R

Question

В настоящее время я пытаюсь выполнить историческую декомпозицию моих рядов данных в R.

Я прочитал тонну статей, и все они дают следующее объяснение того, как сделать историческую декомпозицию:

Где сумма справа - это "динамический прогноз" или "базовый прогноз" Yt+k, зависящий от информации, доступной в момент времени t. Сумма в левой части - это разница между фактическим рядом и базовой проекцией из-за инноваций в переменных в периоды от t+1 до t + k

Я очень запутался в базовом прогнозе и не уверен, какие данные используются!

Мои попытки.

У меня есть 6 переменных VAR, с 55 наблюдениями. Я получаю структурную форму модели, используя разложение Холецкого. После этого я использую функцию Phi, чтобы получить структурное представление скользящего среднего SVAR. Затем я сохраняю этот массив "Phi", чтобы использовать его позже.

    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)

Я думаю, что я получил информацию, необходимую для базовой проекции, но я понятия не имею, куда идти дальше.

ЦЕЛЬ Что я хочу сделать, это оценить влияние 1-й переменной в VAR на 6-ю переменную в VAR за весь период выборки.

Любая помощь будет оценена.

Answer 1

Я перевел функцию VARhd из набора инструментов Matlab Чеза-Бьянки в R код. Моя функция совместима с VAR функция от vars пакеты в R,

Оригинальная функция в 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
% [e-mail]


%% 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)

}

В качестве входного аргумента вы должны использовать VAR функция от vars пакеты в R, Функция возвращает 3-мерный массив: количество наблюдений x количество шоков x количество переменных. (Примечание: я не переводил всю функцию, например, пропустил регистр экзогенных переменных.) Чтобы запустить его, вам понадобятся две дополнительные функции, которые также были переведены из Инструментария Бианки:

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:

Answer 2

Историческая декомпозиция действительно обращается к тому, как ошибки одного ряда влияют на другой ряд в VAR. Самый простой способ сделать это - создать массив подобранных ошибок. Отсюда вам понадобится тройной вложенный цикл:

1. Зациклите установленную серию амортизаторов: for (iShock in 1:6)

2. Цикл по времени измерения данного установленного шока, начиная с периода после базового периода: for (iShockPeriod in 1:55)

3. Смоделируйте эффект индивидуальной реализации этого значения шока для остальной части образца: for (iResponsePeriod in iShockPeriod:55)

В результате получается массив 4D с размерами (например) 6x6x55x55. (i,j,k,l) Элементом будет что-то вроде влияния шока для i-го ряда в k-м периоде на j-й ряд в l-м периоде. Когда я писал ранее реализации, обычно имеет смысл суммировать эти измерения по мере необходимости, чтобы не иметь таких больших массивов.

У меня, к сожалению, нет реализации в R чтобы поделиться, но я работаю над одним в Stata. Я обновлю это ссылкой, если скоро получу в презентабельном состоянии.