Как использовать лассо с пакетом Vars в R
Дорогие коллеги-программисты!
Я пытаюсь проанализировать многомерный набор данных (31 переменная, 1100 наблюдений) с помощью векторной авторегрессии со штрафными санкциями.
Поскольку я использую методы, представленные Diebold et. al (2019), чтобы построить сеть взаимосвязей с помощью матриц разложения дисперсии. Я хотел бы использовать их пакет в r: https://www.rdocumentation.org/packages/vars/versions/1.5-3/topics/fevd
Однако этот пакет можно использовать только с обычной оценкой VAR. Я хотел бы использовать штрафную регрессию, такую как LASSO. Итак, как я могу использовать их пакет в R со штрафным VAR?
Что я пробовал? Это пакет Lassovars на github, однако я не могу использовать его в функции fevd(). Он говорит: использует только оценку из класса Vars.
Надеюсь услышать вас снова!
С уважением,
Барт
1 ответ
Предлагаем возможное решение этой очень интересной проблемы на примере данных ниже. Чтобы соответствовать VAR со штрафом, я использовал недавно выпущенный пакет BigVAR. На данный момент у него нет обычных дополнительных функций, таких как создание FEVD, исторической декомпозиции, прогнозных фан-диаграмм и т. Д., Но все необходимые выходные данные из модели сокращенной формы могут быть получены с помощьюcv.BigVAR. Это то, что я сделал в качестве первого шага ниже. Я предоставлю вам возможность тонкой настройки оценки в сокращенной форме, используя функциональные возможности пакета.
# BigVAR ----
library(BigVAR)
library(expm)
library(data.table)
# Create model
data(Y)
p = 4 # lags
k = ncol(Y) # number of equations
N = nrow(Y) - p # number of obs
# Create a Basic VAR-L (Lasso Penalty) with maximum lag order p=4, 10 gridpoints with lambda optimized according to rolling validation of 1-step ahead MSFE
mod1 = constructModel(Y,p,"Basic",gran=c(150,10),RVAR=FALSE,h=1,cv="Rolling",MN=FALSE,verbose=FALSE,IC=TRUE)
results = cv.BigVAR(mod1)
# Get model estimates
A = results@betaPred[,2:ncol(results@betaPred)] # (k x (k*p)) = (3 x (3*4)) coefficient matrix (reduced form)
Sigma = crossprod(resid(varres))/(N-(k*p)-1)
Оттуда мы можем использовать стандартные формулы для вычисления FEVD. Для доступного введения я настоятельно рекомендую главу 4 "Структурный векторный авторегрессионный анализ" Лутца Килиана и Гельмута Люткеполя. Ниже реализовано все, что нам нужно, от IRF в сокращенной форме до MSPE и, наконец, наша формула для вычисления FEVD в R. Это может быть много, и у меня не было времени отполировать этот код, так что это вполне может быть сделано более эффективно, но, надеюсь, все же информативно. Обратите внимание, что я использую стандартное разложение Холецкого для определения VAR, поэтому применяются обычные следствия порядка переменных.
# A number of helper functions ----
# Compute reduced-form IRFs:
compute_Phi = function(p, k, A_comp, n.ahead) {
J = matrix(0,nrow=k,ncol=k*p)
diag(J) = 1
Phi = lapply(1:n.ahead, function(i) {
J %*% (A_comp %^% (i-1)) %*% t(J)
})
return(Phi)
}
# Compute orthogonalized IRFs:
compute_theta_kj_sq = function(Theta, n.ahead) {
theta_kj_sq = lapply(1:n.ahead, function(h) { # loop over h time periods
out = sapply(1:ncol(Theta[[h]]), function(k) {
terms_to_sum = lapply(1:h, function(i) {
Theta[[i]][k,]**2
})
theta_kj_sq_h = Reduce(`+`, terms_to_sum)
})
colnames(out) = colnames(Theta[[h]])
return(out)
})
return(theta_kj_sq)
}
# Compute mean squared prediction error:
compute_mspe = function(Theta, n.ahead=10) {
mspe = lapply(1:n.ahead, function(h) {
terms_to_sum = lapply(1:h, function(i) {
tcrossprod(Theta[[i]])
})
mspe_h = Reduce(`+`, terms_to_sum)
})
}
# Function for FEVD:
fevd = function(A, Sigma, n.ahead) {
k = dim(A)[1] # number of equations
p = dim(A)[2]/k # number of lags
# Turn into companion form:
if (p>1) {
A_comp = VarptoVar1MC(A,p,k)
} else {
A_comp = A
}
# Compute MSPE: ----
Phi = compute_Phi(p,k,A_comp,n.ahead)
P = t(chol.default(Sigma)) # lower triangular Cholesky factor
B_0 = solve(P) # identification matrix
Theta = lapply(1:length(Phi), function(i) {
Phi[[i]] %*% solve(B_0)
})
theta_kj_sq = compute_theta_kj_sq(Theta, n.ahead) # Squared orthogonaliyed IRFs
mspe = compute_mspe(Theta, n.ahead)
# Compute percentage contributions (i.e. FEVDs): ----
fevd_list = lapply(1:k, function(k) {
t(sapply(1:length(mspe), function(h) {
mspe_k = mspe[[h]][k,k]
theta_k_sq = theta_kj_sq[[h]][,k]
fevd = theta_k_sq/mspe_k
}))
})
# Tidy up
fevd_tidy = data.table::rbindlist(
lapply(1:length(fevd_list), function(k) {
fevd_k = data.table::melt(data.table::data.table(fevd_list[[k]])[,h:=.I], id.vars = "h", variable.name = "j")
fevd_k[,k:=paste0("V",k)]
data.table::setcolorder(fevd_k, c("k", "j", "h"))
})
)
return(fevd_tidy)
}
Наконец, давайте реализуем формулу для n.ahead=20 и нанесите на график результаты.
fevd_res = fevd(A, Sigma, 20)
library(ggplot2)
p = ggplot(data=fevd_res) +
geom_area(aes(x=h, y=value, fill=j)) +
facet_wrap(k ~ .) +
scale_x_continuous(
expand=c(0,0)
) +
scale_fill_discrete(
name = "Shock"
) +
labs(
y = "Variance contribution",
x = "Forecast horizon"
) +
theme_bw()
p
Надеюсь, это будет полезно, и, пожалуйста, кричите, если возникнут дополнительные вопросы.
Наконец, предостережение: я протестировал функцию FEVD, сравнив ее с результатами для стандартного VAR с использованием пакета vars, который вы упомянули в своем вопросе, и он прошел проверку (см. Ниже). Но на этом мое "модульное тестирование" не прошло. Код еще никем не проверялся, так что просто помните об этом и, возможно, самостоятельно исследуйте формулы. Если вы или кто-либо другой заметите какие-либо ошибки, я буду благодарен за отзыв.
Редактировать 1
Для полноты картины добавляем быстрое сравнение результатов, возвращаемых vars::fevd и настраиваемая функция сверху.
# Compare to vars package ----
library(vars)
p = 4 # lags
k = ncol(Y) # number of equations
N = nrow(Y) - p
colnames(Y) = sprintf("V%i", 1:ncol(Y))
n.ahead = 20
varres = vars::VAR(Y,p) # reduced-form model using package command; vars:: to make clear that pkg
# Get estimates for custom fevd function:
Sigma = crossprod(resid(varres))/(N-(k*p)-1)
A = t(
sapply(coef(varres), function(i) {
i[,1]
})
)
A = A[,1:(ncol(A)-1)]
# Run the two different functions:
fevd_pkg = vars::fevd(varres, n.ahead)
fevd_cus = fevd(A, Sigma, n.ahead)
Теперь сравним результаты для первой переменной:
> # Package:
> head(fevd_pkg$V1)
V1 V2 V3
[1,] 1.0000000 0.00000000 0.00000000
[2,] 0.9399842 0.01303013 0.04698572
[3,] 0.9422918 0.01062750 0.04708065
[4,] 0.9231440 0.01409313 0.06276291
[5,] 0.9305901 0.01335727 0.05605267
[6,] 0.9093144 0.01278727 0.07789833
> # Custom:
> head(dcast(fevd_cus[k=="V1"], k+h~j, value.var = "value"))
k h V1 V2 V3
1: V1 1 1.0000000 0.00000000 0.00000000
2: V1 2 0.9399842 0.01303013 0.04698572
3: V1 3 0.9422918 0.01062750 0.04708065
4: V1 4 0.9231440 0.01409313 0.06276291
5: V1 5 0.9305901 0.01335727 0.05605267
6: V1 6 0.9093144 0.01278727 0.07789833
Редактировать 2
Чтобы получить обобщенную ОФВД в R вне зависимости от выходных данных vars::VAR(), вы можете использовать функцию ниже. Я переработал часть исходного кодаfrequencyConnectedness::genFEVD.
# Generalized FEVD ----
library(frequencyConnectedness)
genFEVD_cus = function(
A,
Sigma,
n.ahead,
no.corr=F
) {
k = dim(A)[1] # number of equations
p = dim(A)[2]/k # number of lags
# Turn into companion form:
if (p>1) {
A_comp = BigVAR::VarptoVar1MC(A,p,k)
} else {
A_comp = A
}
# Set off-diagonals to zero:
if (no.corr) {
Sigma = diag(diag(Sigma))
}
Phi = compute_Phi(p,k,A_comp,n.ahead+1) # Reduced-form irfs
denom = diag(Reduce("+", lapply(Phi, function(i) i %*% Sigma %*%
t(i))))
enum = Reduce("+", lapply(Phi, function(i) (i %*% Sigma)^2))
tab = sapply(1:nrow(enum), function(j) enum[j, ]/(denom[j] *
diag(Sigma)))
tab = t(apply(tab, 2, function(i) i/sum(i)))
return(tab)
}
Продолжая предыдущий пример (Правка 1), сравните результаты пользовательской функции с результатамиfrequencyConnectedness::genFEVD который зависит от вывода из vars::VAR():
> frequencyConnectedness::genFEVD(varres, n.ahead)
V1 V2 V3
[1,] 0.89215734 0.02281892 0.08502374
[2,] 0.72025050 0.19335481 0.08639469
[3,] 0.08328267 0.11769438 0.79902295
> genFEVD_cus(A, Sigma, n.ahead)
V1 V2 V3
[1,] 0.89215734 0.02281892 0.08502374
[2,] 0.72025050 0.19335481 0.08639469
[3,] 0.08328267 0.11769438 0.79902295