Как рассчитать доверительные интервалы для этого импульсного отклика? - PullRequest
Новая
       36

Как рассчитать доверительные интервалы для этого импульсного отклика?

1 голос
/ 28 февраля 2020

Я пытаюсь создать доверительные интервалы для моих импульсных откликов. По сути, мне нужен доверительный интервал для каждого наблюдения импульсной характеристики каждой переменной. Мне проще показать вам пример, где я застрял. 

library(vars)
library(varexternalinstrument)
data(GKdata)

gkvar <- VAR(GKdata[, c("logip", "logcpi", "gs1", "ebp")], p = 12, type = "const")
shockcol <- externalinstrument(gkvar, GKdata$ff4_tc, "gs1")
shockcol

ma_representation <- Phi(gkvar, 50)
irfs <- apply(ma_representation, 3, function(x) x %*% shockcol) #unclear
irfs <- as.data.frame(t(irfs))
colnames(irfs) <- names(shockcol)
irfs

Пока все хорошо, мы получаем импульсные отклики для 4 переменных в наборе данных. Теперь для каждого наблюдения в каждой переменной я хочу вычислить начальные доверительные интервалы. Тем не менее, я застрял, вероятно, на статистике для вставки в функцию. Я указал вопросительным знаком под отсутствующим кодом, что я не могу написать 

? <- function() # I need to create a function that says bootstrap for each observation of each variable, I don't know how to do it

boot <- boot(irfs, ?, R = 1000, stype = "w", sim = "ordinary")

boot_ci <- boot.ci(boot, conf = c(0.90, 0.95), type = c("norm", "basic", "perc", "bca"))

Может ли кто-нибудь помочь мне с этим?

Большое спасибо!

1 Ответ

1 голос
/ 28 февраля 2020

Таким образом, данные собираются в период с 1979: 7 по 2012: 6, следовательно, каждые 12 строк - это год, и bootstrap должен выбирать год в блоке.

Сначала мы добавим это информация к данным: 

Data = GKdata
Data$year = (1:nrow(Data) - 1) %/% 12

И мы создаем функцию, которая берет список с выборочными индексами, объединяет их в data.frame и выполняет ту же функцию подгонки. Для загрузки вы не можете вернуть data.frame, поэтому мы возвращаем матрицу: 

func = function(mylist,ind){

  mydf = do.call(rbind,mylist[ind])
  gkvar <- VAR(mydf[, c("logip", "logcpi", "gs1", "ebp")], p = 12, type = "const")
  shockcol <- externalinstrument(gkvar, mydf$ff4_tc, "gs1")

  ma_representation <- Phi(gkvar, 50)
  irfs <- t(apply(ma_representation, 3, function(x) x %*% shockcol)) #unclear
  colnames(irfs) <- names(shockcol)
  irfs
}

Таким образом, мы делаем bootstrap, передавая ему список, разделенный на год. При начальной загрузке будут взяты индексы из этого списка, объединены выбранные годы в новый фрейм данных и произведена подгонка: 

res = boot(split(Data,Data$year),func,100)

Мы можем проверить, правильно ли функция выполняет исходное наблюдение: 

table(res$t0 == irfs)

У ваших irfs есть значения 51 x 4. Если мы посмотрим на оценки в результатах загрузки, это выровняется: 

dim(res$t)
[1] 100 204

Чтобы получить 95% доверительный интервал для каждого значения или переменной, мы делаем: 

boot_est= sapply(1:ncol(res$t),function(i){
  boot.ci(res,index=i,type = "perc")$percent[4:5]
})

А затем создайте верхнюю и нижнюю границу, поместив их обратно в матрицу: 

lb = matrix(boot_est[1,],ncol=ncol(res$t0))
ub = matrix(boot_est[2,],ncol=ncol(res$t0))

Вы можете построить, скажем, первый столбец irfs, почему-то значение ci действительно велико для первых нескольких записей: 

plot(irfs[,1],ylim=c(-8,8))
lines(ub[,1])
lines(lb[,1])

enter image description here

...