R/stationarise_VAR.R
In nicholasjclark/mvgam: Multivariate (Dynamic) Generalized Additive Models
Defines functions
stationarise_VARcor
stationarise_VAR
#### Modifications to Stan code for stationary VAR1 processes ####
# All functions and reparameterisations use code supplied generously by Sarah Heaps:
# Heaps, Sarah E. "Enforcing stationarity through the prior in vector autoregressions."
# Journal of Computational and Graphical Statistics (2022): 1-10.
#' @noRd
stationarise_VAR = function(model_file) {
# Remove previous prior for VAR coefficients
model_file <- model_file[
-c(
(grep('// VAR coefficients', model_file, fixed = TRUE)):(grep(
'// VAR coefficients',
model_file,
fixed = TRUE
) +
1)
)
]
model_file <- model_file[
-c(
(grep("// latent trend VAR1 terms", model_file, fixed = TRUE)):(grep(
"// latent trend VAR1 terms",
model_file,
fixed = TRUE
) +
1)
)
]
model_file <- model_file[
-grep("matrix[n_series, n_series] Sigma;", model_file, fixed = TRUE)
]
model_file <- model_file[
-grep("Sigma = diag_matrix(square(sigma));", model_file, fixed = TRUE)
]
# Add Heaps' functions for constrained VAR1 process priors
if (any(grepl('functions {', model_file, fixed = TRUE))) {
model_file[grep('functions {', model_file, fixed = TRUE)] <-
paste0(
'functions {\n',
'/* Function to compute the matrix square root */\n',
'/* see Heaps 2022 for details (https://doi.org/10.1080/10618600.2022.2079648)*/\n',
'matrix sqrtm(matrix A) {\n',
'int m = rows(A);\n',
'vector[m] root_root_evals = sqrt(sqrt(eigenvalues_sym(A)));\n',
'matrix[m, m] evecs = eigenvectors_sym(A);\n',
'matrix[m, m] eprod = diag_post_multiply(evecs, root_root_evals);\n',
'return tcrossprod(eprod);\n',
'}\n',
'/* Function to transform P_real to P */\n',
'/* see Heaps 2022 for details (https://doi.org/10.1080/10618600.2022.2079648)*/\n',
'matrix P_realtoP(matrix P_real) {\n',
'int m = rows(P_real);\n',
'matrix[m, m] B = tcrossprod(P_real);\n',
'for(i in 1:m) B[i, i] += 1.0;\n',
'return mdivide_left_spd(sqrtm(B), P_real);\n',
'}\n',
'/* Function to perform the reverse mapping*/\n',
'/* see Heaps 2022 for details (https://doi.org/10.1080/10618600.2022.2079648)*/\n',
'matrix[,] rev_mapping(matrix[] P, matrix Sigma) {\n',
'int p = size(P);\n',
'int m = rows(Sigma);\n',
'matrix[m, m] phi_for[p, p]; matrix[m, m] phi_rev[p, p];\n',
'matrix[m, m] Sigma_for[p+1]; matrix[m, m] Sigma_rev[p+1];\n',
'matrix[m, m] S_for; matrix[m, m] S_rev;\n',
'matrix[m, m] S_for_list[p+1];\n',
'matrix[m, m] Gamma_trans[p+1];\n',
'matrix[m, m] phiGamma[2, p];\n',
'// Step 1:\n',
'Sigma_for[p+1] = Sigma;\n',
'S_for_list[p+1] = sqrtm(Sigma);\n',
'for(s in 1:p) {\n',
'// In this block of code S_rev is B^{-1} and S_for is a working matrix\n',
'S_for = - tcrossprod(P[p-s+1]);\n',
'for(i in 1:m) S_for[i, i] += 1.0;\n',
'S_rev = sqrtm(S_for);\n',
'S_for_list[p-s+1] = mdivide_right_spd(mdivide_left_spd(S_rev,\n',
'sqrtm(quad_form_sym(Sigma_for[p-s+2], S_rev))), S_rev);\n',
'Sigma_for[p-s+1] = tcrossprod(S_for_list[p-s+1]);\n',
'}\n',
'// Step 2:\n',
'Sigma_rev[1] = Sigma_for[1];\n',
'Gamma_trans[1] = Sigma_for[1];\n',
'for(s in 0:(p-1)) {\n',
'S_for = S_for_list[s+1];\n',
'S_rev = sqrtm(Sigma_rev[s+1]);\n',
'phi_for[s+1, s+1] = mdivide_right_spd(S_for * P[s+1], S_rev);\n',
"phi_rev[s+1, s+1] = mdivide_right_spd(S_rev * P[s+1]', S_for);\n",
'Gamma_trans[s+2] = phi_for[s+1, s+1] * Sigma_rev[s+1];\n',
'if(s>=1) {\n',
'for(k in 1:s) {\n',
'phi_for[s+1, k] = phi_for[s, k] - phi_for[s+1, s+1] * phi_rev[s, s-k+1];\n',
'phi_rev[s+1, k] = phi_rev[s, k] - phi_rev[s+1, s+1] * phi_for[s, s-k+1];\n',
'}\n',
'for(k in 1:s) Gamma_trans[s+2] = Gamma_trans[s+2] + phi_for[s, k] * Gamma_trans[s+2-k];\n',
'}\n',
"Sigma_rev[s+2] = Sigma_rev[s+1] - quad_form_sym(Sigma_for[s+1],phi_rev[s+1, s+1]');\n",
'}\n',
'for(i in 1:p) phiGamma[1, i] = phi_for[p, i];\n',
"for(i in 1:p) phiGamma[2, i] = Gamma_trans[i]';\n",
'return phiGamma;\n',
'}\n'
)
} else {
model_file[grep('Stan model code', model_file)] <-
paste0(
'// Stan model code generated by package mvgam\n',
'functions {\n',
'/* Function to compute the matrix square root */\n',
'/* see Heaps 2022 for details (https://doi.org/10.1080/10618600.2022.2079648)*/\n',
'matrix sqrtm(matrix A) {\n',
'int m = rows(A);\n',
'vector[m] root_root_evals = sqrt(sqrt(eigenvalues_sym(A)));\n',
'matrix[m, m] evecs = eigenvectors_sym(A);\n',
'matrix[m, m] eprod = diag_post_multiply(evecs, root_root_evals);\n',
'return tcrossprod(eprod);\n',
'}\n',
'/* Function to transform P_real to P */\n',
'/* see Heaps 2022 for details (https://doi.org/10.1080/10618600.2022.2079648)*/\n',
'matrix P_realtoP(matrix P_real) {\n',
'int m = rows(P_real);\n',
'matrix[m, m] B = tcrossprod(P_real);\n',
'for(i in 1:m) B[i, i] += 1.0;\n',
'return mdivide_left_spd(sqrtm(B), P_real);\n',
'}\n',
'/* Function to perform the reverse mapping*/\n',
'/* see Heaps 2022 for details (https://doi.org/10.1080/10618600.2022.2079648)*/\n',
'matrix[,] rev_mapping(matrix[] P, matrix Sigma) {\n',
'int p = size(P);\n',
'int m = rows(Sigma);\n',
'matrix[m, m] phi_for[p, p]; matrix[m, m] phi_rev[p, p];\n',
'matrix[m, m] Sigma_for[p+1]; matrix[m, m] Sigma_rev[p+1];\n',
'matrix[m, m] S_for; matrix[m, m] S_rev;\n',
'matrix[m, m] S_for_list[p+1];\n',
'matrix[m, m] Gamma_trans[p+1];\n',
'matrix[m, m] phiGamma[2, p];\n',
'// Step 1:\n',
'Sigma_for[p+1] = Sigma;\n',
'S_for_list[p+1] = sqrtm(Sigma);\n',
'for(s in 1:p) {\n',
'// In this block of code S_rev is B^{-1} and S_for is a working matrix\n',
'S_for = - tcrossprod(P[p-s+1]);\n',
'for(i in 1:m) S_for[i, i] += 1.0;\n',
'S_rev = sqrtm(S_for);\n',
'S_for_list[p-s+1] = mdivide_right_spd(mdivide_left_spd(S_rev,\n',
'sqrtm(quad_form_sym(Sigma_for[p-s+2], S_rev))), S_rev);\n',
'Sigma_for[p-s+1] = tcrossprod(S_for_list[p-s+1]);\n',
'}\n',
'// Step 2:\n',
'Sigma_rev[1] = Sigma_for[1];\n',
'Gamma_trans[1] = Sigma_for[1];\n',
'for(s in 0:(p-1)) {\n',
'S_for = S_for_list[s+1];\n',
'S_rev = sqrtm(Sigma_rev[s+1]);\n',
'phi_for[s+1, s+1] = mdivide_right_spd(S_for * P[s+1], S_rev);\n',
"phi_rev[s+1, s+1] = mdivide_right_spd(S_rev * P[s+1]', S_for);\n",
'Gamma_trans[s+2] = phi_for[s+1, s+1] * Sigma_rev[s+1];\n',
'if(s>=1) {\n',
'for(k in 1:s) {\n',
'phi_for[s+1, k] = phi_for[s, k] - phi_for[s+1, s+1] * phi_rev[s, s-k+1];\n',
'phi_rev[s+1, k] = phi_rev[s, k] - phi_rev[s+1, s+1] * phi_for[s, s-k+1];\n',
'}\n',
'for(k in 1:s) Gamma_trans[s+2] = Gamma_trans[s+2] + phi_for[s, k] * Gamma_trans[s+2-k];\n',
'}\n',
"Sigma_rev[s+2] = Sigma_rev[s+1] - quad_form_sym(Sigma_for[s+1],phi_rev[s+1, s+1]');\n",
'}\n',
'for(i in 1:p) phiGamma[1, i] = phi_for[p, i];\n',
"for(i in 1:p) phiGamma[2, i] = Gamma_trans[i]';\n",
'return phiGamma;\n',
'}\n',
'}'
)
}
model_file <- readLines(textConnection(model_file), n = -1)
# Add transformed data lines for Heaps stationarity constraints
if (any(grepl('transformed data {', model_file, fixed = TRUE))) {
model_file[grep('transformed data {', model_file, fixed = TRUE)] <- paste0(
'transformed data {\n',
'vector[n_series] trend_zeros = rep_vector(0.0, n_series);\n',
'// exchangeable partial autocorrelation hyperparameters\n',
'// see Heaps 2022 for details (https://doi.org/10.1080/10618600.2022.2079648)\n',
'vector[2] es;\n',
'vector<lower=0>[2] fs;\n',
'vector<lower=0>[2] gs;\n',
'vector<lower=0>[2] hs;\n',
'es[1] = 0;\nes[2] = 0;\n',
'fs[1] = sqrt(0.455);\nfs[2] = sqrt(0.455);\n',
'gs[1] = 1.365;\ngs[2] = 1.365;\n',
'hs[1] = 0.071175;\nhs[2] = 0.071175;\n'
)
} else {
params_line <- min(which(grepl('parameters {', model_file, fixed = TRUE)))
model_file[params_line] <- paste0(
'transformed data {\n',
'vector[n_series] trend_zeros = rep_vector(0.0, n_series);\n',
'// exchangeable partial autocorrelation hyperparameters\n',
'// see Heaps 2022 for details (https://doi.org/10.1080/10618600.2022.2079648)\n',
'vector[2] es;\n',
'vector<lower=0>[2] fs;\n',
'vector<lower=0>[2] gs;\n',
'vector<lower=0>[2] hs;\n',
'es[1] = 0;\nes[2] = 0;\n',
'fs[1] = sqrt(0.455);\nfs[2] = sqrt(0.455);\n',
'gs[1] = 1.365;\ngs[2] = 1.365;\n',
'hs[1] = 0.071175;\nhs[2] = 0.071175;\n',
'}\n\nparameters {'
)
model_file <- readLines(textConnection(model_file), n = -1)
}
# Add parameters for real-valued partial autocorrelations
model_file[grep(
"vector<lower=0>[n_series] sigma;",
model_file,
fixed = TRUE
)] <- paste0(
"vector<lower=0>[n_series] sigma;\n\n",
'// unconstrained VAR1 partial autocorrelations\n',
'matrix[n_series, n_series] P_real;\n',
'// partial autocorrelation hyperparameters\n',
'vector[2] Pmu;\n',
'vector<lower=0>[2] Pomega;\n'
)
model_file <- readLines(textConnection(model_file), n = -1)
# Add transformed parameters for partial autocorrelation reverse mapping
model_file[grep(
'transformed parameters {',
model_file,
fixed = TRUE
)] <- paste0(
'transformed parameters {\n',
"// latent trend VAR1 autoregressive terms\n",
"matrix[n_series, n_series] A;\n",
'// stationary trend covariance\n',
"matrix[n_series, n_series] Sigma;"
)
model_file[
grep(
'trend[i, 1:n_series] = to_row_vector(trend_raw[i]);',
model_file,
fixed = TRUE
) +
1
] <- paste0(
'}\n\n',
"Sigma = diag_matrix(square(sigma));\n",
'// stationary VAR reparameterisation\n',
'{\n',
'matrix[n_series, n_series] P[1];\n',
'matrix[n_series, n_series] phiGamma[2, 1];\n',
'P[1] = P_realtoP(P_real);\n',
'phiGamma = rev_mapping(P, Sigma);\n',
'A = phiGamma[1, 1];\n',
'}'
)
model_file <- readLines(textConnection(model_file), n = -1)
# Add priors for partial autocorrelations
model_file[grep("// trend means", model_file, fixed = TRUE)] <- paste0(
'// partial autocorrelation hyperpriors\n',
'Pmu ~ normal(es, fs);\n',
'Pomega ~ gamma(gs, hs);\n',
'// unconstrained partial autocorrelations\n',
'diagonal(P_real) ~ normal(Pmu[1], 1 / sqrt(Pomega[1]));\n',
'for(i in 1:n_series) {\n',
'for(j in 1:n_series) {\n',
'if(i != j) P_real[i, j] ~ normal(Pmu[2], 1 / sqrt(Pomega[2]));\n',
'}\n',
'}\n',
"// trend means"
)
model_file <- readLines(textConnection(model_file), n = -1)
# Return the updated model file
return(model_file)
}
#' Modifications for a VAR1 with possible contemporaneously correlated trend
#' errors
#' @noRd
stationarise_VARcor = function(model_file) {
# Remove previous prior for VAR coefficients
model_file <- model_file[
-c(
(grep('// VAR coefficients', model_file, fixed = TRUE)):(grep(
'// VAR coefficients',
model_file,
fixed = TRUE
) +
1)
)
]
model_file <- model_file[
-c(
(grep("// latent trend VAR1 terms", model_file, fixed = TRUE)):(grep(
"// latent trend VAR1 terms",
model_file,
fixed = TRUE
) +
1)
)
]
model_file <- model_file[
-grep("matrix[n_series, n_series] Sigma;", model_file, fixed = TRUE)
]
model_file <- model_file[
-grep("Sigma = diag_matrix(square(sigma));", model_file, fixed = TRUE)
]
# Add Heaps' functions for constrained VAR1 process priors
if (any(grepl('functions {', model_file, fixed = TRUE))) {
model_file[grep('functions {', model_file, fixed = TRUE)] <-
paste0(
'functions {\n',
'/* Function to compute the matrix square root */\n',
'/* see Heaps 2022 for details (https://doi.org/10.1080/10618600.2022.2079648)*/\n',
'matrix sqrtm(matrix A) {\n',
'int m = rows(A);\n',
'vector[m] root_root_evals = sqrt(sqrt(eigenvalues_sym(A)));\n',
'matrix[m, m] evecs = eigenvectors_sym(A);\n',
'matrix[m, m] eprod = diag_post_multiply(evecs, root_root_evals);\n',
'return tcrossprod(eprod);\n',
'}\n',
'/* Function to transform P_real to P */\n',
'/* see Heaps 2022 for details (https://doi.org/10.1080/10618600.2022.2079648)*/\n',
'matrix P_realtoP(matrix P_real) {\n',
'int m = rows(P_real);\n',
'matrix[m, m] B = tcrossprod(P_real);\n',
'for(i in 1:m) B[i, i] += 1.0;\n',
'return mdivide_left_spd(sqrtm(B), P_real);\n',
'}\n',
'/* Function to perform the reverse mapping*/\n',
'/* see Heaps 2022 for details (https://doi.org/10.1080/10618600.2022.2079648)*/\n',
'matrix[,] rev_mapping(matrix[] P, matrix Sigma) {\n',
'int p = size(P);\n',
'int m = rows(Sigma);\n',
'matrix[m, m] phi_for[p, p]; matrix[m, m] phi_rev[p, p];\n',
'matrix[m, m] Sigma_for[p+1]; matrix[m, m] Sigma_rev[p+1];\n',
'matrix[m, m] S_for; matrix[m, m] S_rev;\n',
'matrix[m, m] S_for_list[p+1];\n',
'matrix[m, m] Gamma_trans[p+1];\n',
'matrix[m, m] phiGamma[2, p];\n',
'// Step 1:\n',
'Sigma_for[p+1] = Sigma;\n',
'S_for_list[p+1] = sqrtm(Sigma);\n',
'for(s in 1:p) {\n',
'// In this block of code S_rev is B^{-1} and S_for is a working matrix\n',
'S_for = - tcrossprod(P[p-s+1]);\n',
'for(i in 1:m) S_for[i, i] += 1.0;\n',
'S_rev = sqrtm(S_for);\n',
'S_for_list[p-s+1] = mdivide_right_spd(mdivide_left_spd(S_rev,\n',
'sqrtm(quad_form_sym(Sigma_for[p-s+2], S_rev))), S_rev);\n',
'Sigma_for[p-s+1] = tcrossprod(S_for_list[p-s+1]);\n',
'}\n',
'// Step 2:\n',
'Sigma_rev[1] = Sigma_for[1];\n',
'Gamma_trans[1] = Sigma_for[1];\n',
'for(s in 0:(p-1)) {\n',
'S_for = S_for_list[s+1];\n',
'S_rev = sqrtm(Sigma_rev[s+1]);\n',
'phi_for[s+1, s+1] = mdivide_right_spd(S_for * P[s+1], S_rev);\n',
"phi_rev[s+1, s+1] = mdivide_right_spd(S_rev * P[s+1]', S_for);\n",
'Gamma_trans[s+2] = phi_for[s+1, s+1] * Sigma_rev[s+1];\n',
'if(s>=1) {\n',
'for(k in 1:s) {\n',
'phi_for[s+1, k] = phi_for[s, k] - phi_for[s+1, s+1] * phi_rev[s, s-k+1];\n',
'phi_rev[s+1, k] = phi_rev[s, k] - phi_rev[s+1, s+1] * phi_for[s, s-k+1];\n',
'}\n',
'for(k in 1:s) Gamma_trans[s+2] = Gamma_trans[s+2] + phi_for[s, k] * Gamma_trans[s+2-k];\n',
'}\n',
"Sigma_rev[s+2] = Sigma_rev[s+1] - quad_form_sym(Sigma_for[s+1],phi_rev[s+1, s+1]');\n",
'}\n',
'for(i in 1:p) phiGamma[1, i] = phi_for[p, i];\n',
"for(i in 1:p) phiGamma[2, i] = Gamma_trans[i]';\n",
'return phiGamma;\n',
'}\n'
)
} else {
model_file[grep('Stan model code', model_file)] <-
paste0(
'// Stan model code generated by package mvgam\n',
'functions {\n',
'/* Function to compute the matrix square root */\n',
'/* see Heaps 2022 for details (https://doi.org/10.1080/10618600.2022.2079648)*/\n',
'matrix sqrtm(matrix A) {\n',
'int m = rows(A);\n',
'vector[m] root_root_evals = sqrt(sqrt(eigenvalues_sym(A)));\n',
'matrix[m, m] evecs = eigenvectors_sym(A);\n',
'matrix[m, m] eprod = diag_post_multiply(evecs, root_root_evals);\n',
'return tcrossprod(eprod);\n',
'}\n',
'/* Function to transform P_real to P */\n',
'/* see Heaps 2022 for details (https://doi.org/10.1080/10618600.2022.2079648)*/\n',
'matrix P_realtoP(matrix P_real) {\n',
'int m = rows(P_real);\n',
'matrix[m, m] B = tcrossprod(P_real);\n',
'for(i in 1:m) B[i, i] += 1.0;\n',
'return mdivide_left_spd(sqrtm(B), P_real);\n',
'}\n',
'/* Function to perform the reverse mapping*/\n',
'/* see Heaps 2022 for details (https://doi.org/10.1080/10618600.2022.2079648)*/\n',
'matrix[,] rev_mapping(matrix[] P, matrix Sigma) {\n',
'int p = size(P);\n',
'int m = rows(Sigma);\n',
'matrix[m, m] phi_for[p, p]; matrix[m, m] phi_rev[p, p];\n',
'matrix[m, m] Sigma_for[p+1]; matrix[m, m] Sigma_rev[p+1];\n',
'matrix[m, m] S_for; matrix[m, m] S_rev;\n',
'matrix[m, m] S_for_list[p+1];\n',
'matrix[m, m] Gamma_trans[p+1];\n',
'matrix[m, m] phiGamma[2, p];\n',
'// Step 1:\n',
'Sigma_for[p+1] = Sigma;\n',
'S_for_list[p+1] = sqrtm(Sigma);\n',
'for(s in 1:p) {\n',
'// In this block of code S_rev is B^{-1} and S_for is a working matrix\n',
'S_for = - tcrossprod(P[p-s+1]);\n',
'for(i in 1:m) S_for[i, i] += 1.0;\n',
'S_rev = sqrtm(S_for);\n',
'S_for_list[p-s+1] = mdivide_right_spd(mdivide_left_spd(S_rev,\n',
'sqrtm(quad_form_sym(Sigma_for[p-s+2], S_rev))), S_rev);\n',
'Sigma_for[p-s+1] = tcrossprod(S_for_list[p-s+1]);\n',
'}\n',
'// Step 2:\n',
'Sigma_rev[1] = Sigma_for[1];\n',
'Gamma_trans[1] = Sigma_for[1];\n',
'for(s in 0:(p-1)) {\n',
'S_for = S_for_list[s+1];\n',
'S_rev = sqrtm(Sigma_rev[s+1]);\n',
'phi_for[s+1, s+1] = mdivide_right_spd(S_for * P[s+1], S_rev);\n',
"phi_rev[s+1, s+1] = mdivide_right_spd(S_rev * P[s+1]', S_for);\n",
'Gamma_trans[s+2] = phi_for[s+1, s+1] * Sigma_rev[s+1];\n',
'if(s>=1) {\n',
'for(k in 1:s) {\n',
'phi_for[s+1, k] = phi_for[s, k] - phi_for[s+1, s+1] * phi_rev[s, s-k+1];\n',
'phi_rev[s+1, k] = phi_rev[s, k] - phi_rev[s+1, s+1] * phi_for[s, s-k+1];\n',
'}\n',
'for(k in 1:s) Gamma_trans[s+2] = Gamma_trans[s+2] + phi_for[s, k] * Gamma_trans[s+2-k];\n',
'}\n',
"Sigma_rev[s+2] = Sigma_rev[s+1] - quad_form_sym(Sigma_for[s+1],phi_rev[s+1, s+1]');\n",
'}\n',
'for(i in 1:p) phiGamma[1, i] = phi_for[p, i];\n',
"for(i in 1:p) phiGamma[2, i] = Gamma_trans[i]';\n",
'return phiGamma;\n',
'}\n',
'}'
)
}
model_file <- readLines(textConnection(model_file), n = -1)
# Add transformed data lines for Heaps stationarity constraints
if (any(grepl('transformed data {', model_file, fixed = TRUE))) {
model_file[grep('transformed data {', model_file, fixed = TRUE)] <- paste0(
'transformed data {\n',
'vector[n_series] trend_zeros = rep_vector(0.0, n_series);\n',
'// exchangeable partial autocorrelation hyperparameters\n',
'// see Heaps 2022 for details (https://doi.org/10.1080/10618600.2022.2079648)\n',
'vector[2] es;\n',
'vector<lower=0>[2] fs;\n',
'vector<lower=0>[2] gs;\n',
'vector<lower=0>[2] hs;\n',
'es[1] = 0;\nes[2] = 0;\n',
'fs[1] = sqrt(0.455);\nfs[2] = sqrt(0.455);\n',
'gs[1] = 1.365;\ngs[2] = 1.365;\n',
'hs[1] = 0.071175;\nhs[2] = 0.071175;\n'
)
} else {
params_line <- min(which(grepl('parameters {', model_file, fixed = TRUE)))
model_file[params_line] <- paste0(
'transformed data {\n',
'vector[n_series] trend_zeros = rep_vector(0.0, n_series);\n',
'// exchangeable partial autocorrelation hyperparameters\n',
'// see Heaps 2022 for details (https://doi.org/10.1080/10618600.2022.2079648)\n',
'vector[2] es;\n',
'vector<lower=0>[2] fs;\n',
'vector<lower=0>[2] gs;\n',
'vector<lower=0>[2] hs;\n',
'es[1] = 0;\nes[2] = 0;\n',
'fs[1] = sqrt(0.455);\nfs[2] = sqrt(0.455);\n',
'gs[1] = 1.365;\ngs[2] = 1.365;\n',
'hs[1] = 0.071175;\nhs[2] = 0.071175;\n',
'}\n\nparameters {'
)
model_file <- readLines(textConnection(model_file), n = -1)
}
# Add parameters for real-valued partial autocorrelations
model_file[grep(
"vector<lower=0>[n_series] sigma;",
model_file,
fixed = TRUE
)] <- paste0(
"cholesky_factor_corr[n_series] L_Omega;\n",
"vector<lower=0>[n_series] sigma;\n\n",
'// unconstrained VAR1 partial autocorrelations\n',
'matrix[n_series, n_series] P_real;\n',
'// partial autocorrelation hyperparameters\n',
'vector[2] Pmu;\n',
'vector<lower=0>[2] Pomega;\n'
)
model_file <- readLines(textConnection(model_file), n = -1)
# Add transformed parameters for partial autocorrelation reverse mapping
model_file[grep(
'transformed parameters {',
model_file,
fixed = TRUE
)] <- paste0(
'transformed parameters {\n',
"// latent trend VAR1 autoregressive terms\n",
"matrix[n_series, n_series] A;\n",
'// LKJ form of covariance matrix\n',
"matrix[n_series, n_series] L_Sigma;\n",
'// computed error covariance matrix\n',
'cov_matrix[n_series] Sigma;\n',
'// initial trend covariance\n',
'cov_matrix[n_series] Gamma;'
)
model_file[
grep(
'trend[i, 1:n_series] = to_row_vector(trend_raw[i]);',
model_file,
fixed = TRUE
) +
1
] <- paste0(
'}\n\n',
'L_Sigma = diag_pre_multiply(sigma, L_Omega);\n',
'Sigma = multiply_lower_tri_self_transpose(L_Sigma);\n',
'// stationary VAR reparameterisation\n',
'{\n',
'matrix[n_series, n_series] P[1];\n',
'matrix[n_series, n_series] phiGamma[2, 1];\n',
'P[1] = P_realtoP(P_real);\n',
'phiGamma = rev_mapping(P, Sigma);\n',
'A = phiGamma[1, 1];\n',
'Gamma = phiGamma[2, 1];\n',
'}'
)
model_file <- readLines(textConnection(model_file), n = -1)
# Add priors for partial autocorrelations
model_file[grep("// trend means", model_file, fixed = TRUE)] <- paste0(
'// LKJ error correlation prior\n',
'L_Omega ~ lkj_corr_cholesky(2);\n',
'// partial autocorrelation hyperpriors\n',
'Pmu ~ normal(es, fs);\n',
'Pomega ~ gamma(gs, hs);\n',
'// unconstrained partial autocorrelations\n',
'diagonal(P_real) ~ normal(Pmu[1], 1 / sqrt(Pomega[1]));\n',
'for(i in 1:n_series) {\n',
'for(j in 1:n_series) {\n',
'if(i != j) P_real[i, j] ~ normal(Pmu[2], 1 / sqrt(Pomega[2]));\n',
'}\n',
'}\n',
"// trend means"
)
model_file <- readLines(textConnection(model_file), n = -1)
# Update prior for error SD parameters sigma
model_file[
grep(
"// priors for latent trend variance parameters",
model_file,
fixed = TRUE
) +
1
] <- paste0(
'sigma ~ inv_gamma(2.3693353, 0.7311319);'
)
model_file <- readLines(textConnection(model_file), n = -1)
# Update trend model to use multinormal
model_file[grep(
"trend_raw[1] ~ normal(0, sigma);",
model_file,
fixed = TRUE
)] <- "trend_raw[1] ~ multi_normal(trend_zeros, Gamma);"
model_file[grep(
"trend_raw[i] ~ normal(mu[i - 1], sigma);",
model_file,
fixed = TRUE
)] <- "trend_raw[i] ~ multi_normal_cholesky(mu[i - 1], L_Sigma);"
model_file <- gsub(
'contemporaneously uncorrelated',
'contemporaneously correlated',
model_file
)
model_file <- readLines(textConnection(model_file), n = -1)
# Return the updated model file
return(model_file)
}
nicholasjclark/mvgam documentation built on April 17, 2025, 9:39 p.m.
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Defines functions stationarise_VARcor stationarise_VAR
#### Modifications to Stan code for stationary VAR1 processes ####
# All functions and reparameterisations use code supplied generously by Sarah Heaps:
# Heaps, Sarah E. "Enforcing stationarity through the prior in vector autoregressions."
# Journal of Computational and Graphical Statistics (2022): 1-10.
#' @noRd
stationarise_VAR = function(model_file) {
# Remove previous prior for VAR coefficients
model_file <- model_file[
-c(
(grep('// VAR coefficients', model_file, fixed = TRUE)):(grep(
'// VAR coefficients',
model_file,
fixed = TRUE
) +
1)
)
]
model_file <- model_file[
-c(
(grep("// latent trend VAR1 terms", model_file, fixed = TRUE)):(grep(
"// latent trend VAR1 terms",
model_file,
fixed = TRUE
) +
1)
)
]
model_file <- model_file[
-grep("matrix[n_series, n_series] Sigma;", model_file, fixed = TRUE)
]
model_file <- model_file[
-grep("Sigma = diag_matrix(square(sigma));", model_file, fixed = TRUE)
]
# Add Heaps' functions for constrained VAR1 process priors
if (any(grepl('functions {', model_file, fixed = TRUE))) {
model_file[grep('functions {', model_file, fixed = TRUE)] <-
paste0(
'functions {\n',
'/* Function to compute the matrix square root */\n',
'/* see Heaps 2022 for details (https://doi.org/10.1080/10618600.2022.2079648)*/\n',
'matrix sqrtm(matrix A) {\n',
'int m = rows(A);\n',
'vector[m] root_root_evals = sqrt(sqrt(eigenvalues_sym(A)));\n',
'matrix[m, m] evecs = eigenvectors_sym(A);\n',
'matrix[m, m] eprod = diag_post_multiply(evecs, root_root_evals);\n',
'return tcrossprod(eprod);\n',
'}\n',
'/* Function to transform P_real to P */\n',
'/* see Heaps 2022 for details (https://doi.org/10.1080/10618600.2022.2079648)*/\n',
'matrix P_realtoP(matrix P_real) {\n',
'int m = rows(P_real);\n',
'matrix[m, m] B = tcrossprod(P_real);\n',
'for(i in 1:m) B[i, i] += 1.0;\n',
'return mdivide_left_spd(sqrtm(B), P_real);\n',
'}\n',
'/* Function to perform the reverse mapping*/\n',
'/* see Heaps 2022 for details (https://doi.org/10.1080/10618600.2022.2079648)*/\n',
'matrix[,] rev_mapping(matrix[] P, matrix Sigma) {\n',
'int p = size(P);\n',
'int m = rows(Sigma);\n',
'matrix[m, m] phi_for[p, p]; matrix[m, m] phi_rev[p, p];\n',
'matrix[m, m] Sigma_for[p+1]; matrix[m, m] Sigma_rev[p+1];\n',
'matrix[m, m] S_for; matrix[m, m] S_rev;\n',
'matrix[m, m] S_for_list[p+1];\n',
'matrix[m, m] Gamma_trans[p+1];\n',
'matrix[m, m] phiGamma[2, p];\n',
'// Step 1:\n',
'Sigma_for[p+1] = Sigma;\n',
'S_for_list[p+1] = sqrtm(Sigma);\n',
'for(s in 1:p) {\n',
'// In this block of code S_rev is B^{-1} and S_for is a working matrix\n',
'S_for = - tcrossprod(P[p-s+1]);\n',
'for(i in 1:m) S_for[i, i] += 1.0;\n',
'S_rev = sqrtm(S_for);\n',
'S_for_list[p-s+1] = mdivide_right_spd(mdivide_left_spd(S_rev,\n',
'sqrtm(quad_form_sym(Sigma_for[p-s+2], S_rev))), S_rev);\n',
'Sigma_for[p-s+1] = tcrossprod(S_for_list[p-s+1]);\n',
'}\n',
'// Step 2:\n',
'Sigma_rev[1] = Sigma_for[1];\n',
'Gamma_trans[1] = Sigma_for[1];\n',
'for(s in 0:(p-1)) {\n',
'S_for = S_for_list[s+1];\n',
'S_rev = sqrtm(Sigma_rev[s+1]);\n',
'phi_for[s+1, s+1] = mdivide_right_spd(S_for * P[s+1], S_rev);\n',
"phi_rev[s+1, s+1] = mdivide_right_spd(S_rev * P[s+1]', S_for);\n",
'Gamma_trans[s+2] = phi_for[s+1, s+1] * Sigma_rev[s+1];\n',
'if(s>=1) {\n',
'for(k in 1:s) {\n',
'phi_for[s+1, k] = phi_for[s, k] - phi_for[s+1, s+1] * phi_rev[s, s-k+1];\n',
'phi_rev[s+1, k] = phi_rev[s, k] - phi_rev[s+1, s+1] * phi_for[s, s-k+1];\n',
'}\n',
'for(k in 1:s) Gamma_trans[s+2] = Gamma_trans[s+2] + phi_for[s, k] * Gamma_trans[s+2-k];\n',
'}\n',
"Sigma_rev[s+2] = Sigma_rev[s+1] - quad_form_sym(Sigma_for[s+1],phi_rev[s+1, s+1]');\n",
'}\n',
'for(i in 1:p) phiGamma[1, i] = phi_for[p, i];\n',
"for(i in 1:p) phiGamma[2, i] = Gamma_trans[i]';\n",
'return phiGamma;\n',
'}\n'
)
} else {
model_file[grep('Stan model code', model_file)] <-
paste0(
'// Stan model code generated by package mvgam\n',
'functions {\n',
'/* Function to compute the matrix square root */\n',
'/* see Heaps 2022 for details (https://doi.org/10.1080/10618600.2022.2079648)*/\n',
'matrix sqrtm(matrix A) {\n',
'int m = rows(A);\n',
'vector[m] root_root_evals = sqrt(sqrt(eigenvalues_sym(A)));\n',
'matrix[m, m] evecs = eigenvectors_sym(A);\n',
'matrix[m, m] eprod = diag_post_multiply(evecs, root_root_evals);\n',
'return tcrossprod(eprod);\n',
'}\n',
'/* Function to transform P_real to P */\n',
'/* see Heaps 2022 for details (https://doi.org/10.1080/10618600.2022.2079648)*/\n',
'matrix P_realtoP(matrix P_real) {\n',
'int m = rows(P_real);\n',
'matrix[m, m] B = tcrossprod(P_real);\n',
'for(i in 1:m) B[i, i] += 1.0;\n',
'return mdivide_left_spd(sqrtm(B), P_real);\n',
'}\n',
'/* Function to perform the reverse mapping*/\n',
'/* see Heaps 2022 for details (https://doi.org/10.1080/10618600.2022.2079648)*/\n',
'matrix[,] rev_mapping(matrix[] P, matrix Sigma) {\n',
'int p = size(P);\n',
'int m = rows(Sigma);\n',
'matrix[m, m] phi_for[p, p]; matrix[m, m] phi_rev[p, p];\n',
'matrix[m, m] Sigma_for[p+1]; matrix[m, m] Sigma_rev[p+1];\n',
'matrix[m, m] S_for; matrix[m, m] S_rev;\n',
'matrix[m, m] S_for_list[p+1];\n',
'matrix[m, m] Gamma_trans[p+1];\n',
'matrix[m, m] phiGamma[2, p];\n',
'// Step 1:\n',
'Sigma_for[p+1] = Sigma;\n',
'S_for_list[p+1] = sqrtm(Sigma);\n',
'for(s in 1:p) {\n',
'// In this block of code S_rev is B^{-1} and S_for is a working matrix\n',
'S_for = - tcrossprod(P[p-s+1]);\n',
'for(i in 1:m) S_for[i, i] += 1.0;\n',
'S_rev = sqrtm(S_for);\n',
'S_for_list[p-s+1] = mdivide_right_spd(mdivide_left_spd(S_rev,\n',
'sqrtm(quad_form_sym(Sigma_for[p-s+2], S_rev))), S_rev);\n',
'Sigma_for[p-s+1] = tcrossprod(S_for_list[p-s+1]);\n',
'}\n',
'// Step 2:\n',
'Sigma_rev[1] = Sigma_for[1];\n',
'Gamma_trans[1] = Sigma_for[1];\n',
'for(s in 0:(p-1)) {\n',
'S_for = S_for_list[s+1];\n',
'S_rev = sqrtm(Sigma_rev[s+1]);\n',
'phi_for[s+1, s+1] = mdivide_right_spd(S_for * P[s+1], S_rev);\n',
"phi_rev[s+1, s+1] = mdivide_right_spd(S_rev * P[s+1]', S_for);\n",
'Gamma_trans[s+2] = phi_for[s+1, s+1] * Sigma_rev[s+1];\n',
'if(s>=1) {\n',
'for(k in 1:s) {\n',
'phi_for[s+1, k] = phi_for[s, k] - phi_for[s+1, s+1] * phi_rev[s, s-k+1];\n',
'phi_rev[s+1, k] = phi_rev[s, k] - phi_rev[s+1, s+1] * phi_for[s, s-k+1];\n',
'}\n',
'for(k in 1:s) Gamma_trans[s+2] = Gamma_trans[s+2] + phi_for[s, k] * Gamma_trans[s+2-k];\n',
'}\n',
"Sigma_rev[s+2] = Sigma_rev[s+1] - quad_form_sym(Sigma_for[s+1],phi_rev[s+1, s+1]');\n",
'}\n',
'for(i in 1:p) phiGamma[1, i] = phi_for[p, i];\n',
"for(i in 1:p) phiGamma[2, i] = Gamma_trans[i]';\n",
'return phiGamma;\n',
'}\n',
'}'
)
}
model_file <- readLines(textConnection(model_file), n = -1)
# Add transformed data lines for Heaps stationarity constraints
if (any(grepl('transformed data {', model_file, fixed = TRUE))) {
model_file[grep('transformed data {', model_file, fixed = TRUE)] <- paste0(
'transformed data {\n',
'vector[n_series] trend_zeros = rep_vector(0.0, n_series);\n',
'// exchangeable partial autocorrelation hyperparameters\n',
'// see Heaps 2022 for details (https://doi.org/10.1080/10618600.2022.2079648)\n',
'vector[2] es;\n',
'vector<lower=0>[2] fs;\n',
'vector<lower=0>[2] gs;\n',
'vector<lower=0>[2] hs;\n',
'es[1] = 0;\nes[2] = 0;\n',
'fs[1] = sqrt(0.455);\nfs[2] = sqrt(0.455);\n',
'gs[1] = 1.365;\ngs[2] = 1.365;\n',
'hs[1] = 0.071175;\nhs[2] = 0.071175;\n'
)
} else {
params_line <- min(which(grepl('parameters {', model_file, fixed = TRUE)))
model_file[params_line] <- paste0(
'transformed data {\n',
'vector[n_series] trend_zeros = rep_vector(0.0, n_series);\n',
'// exchangeable partial autocorrelation hyperparameters\n',
'// see Heaps 2022 for details (https://doi.org/10.1080/10618600.2022.2079648)\n',
'vector[2] es;\n',
'vector<lower=0>[2] fs;\n',
'vector<lower=0>[2] gs;\n',
'vector<lower=0>[2] hs;\n',
'es[1] = 0;\nes[2] = 0;\n',
'fs[1] = sqrt(0.455);\nfs[2] = sqrt(0.455);\n',
'gs[1] = 1.365;\ngs[2] = 1.365;\n',
'hs[1] = 0.071175;\nhs[2] = 0.071175;\n',
'}\n\nparameters {'
)
model_file <- readLines(textConnection(model_file), n = -1)
}
# Add parameters for real-valued partial autocorrelations
model_file[grep(
"vector<lower=0>[n_series] sigma;",
model_file,
fixed = TRUE
)] <- paste0(
"vector<lower=0>[n_series] sigma;\n\n",
'// unconstrained VAR1 partial autocorrelations\n',
'matrix[n_series, n_series] P_real;\n',
'// partial autocorrelation hyperparameters\n',
'vector[2] Pmu;\n',
'vector<lower=0>[2] Pomega;\n'
)
model_file <- readLines(textConnection(model_file), n = -1)
# Add transformed parameters for partial autocorrelation reverse mapping
model_file[grep(
'transformed parameters {',
model_file,
fixed = TRUE
)] <- paste0(
'transformed parameters {\n',
"// latent trend VAR1 autoregressive terms\n",
"matrix[n_series, n_series] A;\n",
'// stationary trend covariance\n',
"matrix[n_series, n_series] Sigma;"
)
model_file[
grep(
'trend[i, 1:n_series] = to_row_vector(trend_raw[i]);',
model_file,
fixed = TRUE
) +
1
] <- paste0(
'}\n\n',
"Sigma = diag_matrix(square(sigma));\n",
'// stationary VAR reparameterisation\n',
'{\n',
'matrix[n_series, n_series] P[1];\n',
'matrix[n_series, n_series] phiGamma[2, 1];\n',
'P[1] = P_realtoP(P_real);\n',
'phiGamma = rev_mapping(P, Sigma);\n',
'A = phiGamma[1, 1];\n',
'}'
)
model_file <- readLines(textConnection(model_file), n = -1)
# Add priors for partial autocorrelations
model_file[grep("// trend means", model_file, fixed = TRUE)] <- paste0(
'// partial autocorrelation hyperpriors\n',
'Pmu ~ normal(es, fs);\n',
'Pomega ~ gamma(gs, hs);\n',
'// unconstrained partial autocorrelations\n',
'diagonal(P_real) ~ normal(Pmu[1], 1 / sqrt(Pomega[1]));\n',
'for(i in 1:n_series) {\n',
'for(j in 1:n_series) {\n',
'if(i != j) P_real[i, j] ~ normal(Pmu[2], 1 / sqrt(Pomega[2]));\n',
'}\n',
'}\n',
"// trend means"
)
model_file <- readLines(textConnection(model_file), n = -1)
# Return the updated model file
return(model_file)
}
#' Modifications for a VAR1 with possible contemporaneously correlated trend
#' errors
#' @noRd
stationarise_VARcor = function(model_file) {
# Remove previous prior for VAR coefficients
model_file <- model_file[
-c(
(grep('// VAR coefficients', model_file, fixed = TRUE)):(grep(
'// VAR coefficients',
model_file,
fixed = TRUE
) +
1)
)
]
model_file <- model_file[
-c(
(grep("// latent trend VAR1 terms", model_file, fixed = TRUE)):(grep(
"// latent trend VAR1 terms",
model_file,
fixed = TRUE
) +
1)
)
]
model_file <- model_file[
-grep("matrix[n_series, n_series] Sigma;", model_file, fixed = TRUE)
]
model_file <- model_file[
-grep("Sigma = diag_matrix(square(sigma));", model_file, fixed = TRUE)
]
# Add Heaps' functions for constrained VAR1 process priors
if (any(grepl('functions {', model_file, fixed = TRUE))) {
model_file[grep('functions {', model_file, fixed = TRUE)] <-
paste0(
'functions {\n',
'/* Function to compute the matrix square root */\n',
'/* see Heaps 2022 for details (https://doi.org/10.1080/10618600.2022.2079648)*/\n',
'matrix sqrtm(matrix A) {\n',
'int m = rows(A);\n',
'vector[m] root_root_evals = sqrt(sqrt(eigenvalues_sym(A)));\n',
'matrix[m, m] evecs = eigenvectors_sym(A);\n',
'matrix[m, m] eprod = diag_post_multiply(evecs, root_root_evals);\n',
'return tcrossprod(eprod);\n',
'}\n',
'/* Function to transform P_real to P */\n',
'/* see Heaps 2022 for details (https://doi.org/10.1080/10618600.2022.2079648)*/\n',
'matrix P_realtoP(matrix P_real) {\n',
'int m = rows(P_real);\n',
'matrix[m, m] B = tcrossprod(P_real);\n',
'for(i in 1:m) B[i, i] += 1.0;\n',
'return mdivide_left_spd(sqrtm(B), P_real);\n',
'}\n',
'/* Function to perform the reverse mapping*/\n',
'/* see Heaps 2022 for details (https://doi.org/10.1080/10618600.2022.2079648)*/\n',
'matrix[,] rev_mapping(matrix[] P, matrix Sigma) {\n',
'int p = size(P);\n',
'int m = rows(Sigma);\n',
'matrix[m, m] phi_for[p, p]; matrix[m, m] phi_rev[p, p];\n',
'matrix[m, m] Sigma_for[p+1]; matrix[m, m] Sigma_rev[p+1];\n',
'matrix[m, m] S_for; matrix[m, m] S_rev;\n',
'matrix[m, m] S_for_list[p+1];\n',
'matrix[m, m] Gamma_trans[p+1];\n',
'matrix[m, m] phiGamma[2, p];\n',
'// Step 1:\n',
'Sigma_for[p+1] = Sigma;\n',
'S_for_list[p+1] = sqrtm(Sigma);\n',
'for(s in 1:p) {\n',
'// In this block of code S_rev is B^{-1} and S_for is a working matrix\n',
'S_for = - tcrossprod(P[p-s+1]);\n',
'for(i in 1:m) S_for[i, i] += 1.0;\n',
'S_rev = sqrtm(S_for);\n',
'S_for_list[p-s+1] = mdivide_right_spd(mdivide_left_spd(S_rev,\n',
'sqrtm(quad_form_sym(Sigma_for[p-s+2], S_rev))), S_rev);\n',
'Sigma_for[p-s+1] = tcrossprod(S_for_list[p-s+1]);\n',
'}\n',
'// Step 2:\n',
'Sigma_rev[1] = Sigma_for[1];\n',
'Gamma_trans[1] = Sigma_for[1];\n',
'for(s in 0:(p-1)) {\n',
'S_for = S_for_list[s+1];\n',
'S_rev = sqrtm(Sigma_rev[s+1]);\n',
'phi_for[s+1, s+1] = mdivide_right_spd(S_for * P[s+1], S_rev);\n',
"phi_rev[s+1, s+1] = mdivide_right_spd(S_rev * P[s+1]', S_for);\n",
'Gamma_trans[s+2] = phi_for[s+1, s+1] * Sigma_rev[s+1];\n',
'if(s>=1) {\n',
'for(k in 1:s) {\n',
'phi_for[s+1, k] = phi_for[s, k] - phi_for[s+1, s+1] * phi_rev[s, s-k+1];\n',
'phi_rev[s+1, k] = phi_rev[s, k] - phi_rev[s+1, s+1] * phi_for[s, s-k+1];\n',
'}\n',
'for(k in 1:s) Gamma_trans[s+2] = Gamma_trans[s+2] + phi_for[s, k] * Gamma_trans[s+2-k];\n',
'}\n',
"Sigma_rev[s+2] = Sigma_rev[s+1] - quad_form_sym(Sigma_for[s+1],phi_rev[s+1, s+1]');\n",
'}\n',
'for(i in 1:p) phiGamma[1, i] = phi_for[p, i];\n',
"for(i in 1:p) phiGamma[2, i] = Gamma_trans[i]';\n",
'return phiGamma;\n',
'}\n'
)
} else {
model_file[grep('Stan model code', model_file)] <-
paste0(
'// Stan model code generated by package mvgam\n',
'functions {\n',
'/* Function to compute the matrix square root */\n',
'/* see Heaps 2022 for details (https://doi.org/10.1080/10618600.2022.2079648)*/\n',
'matrix sqrtm(matrix A) {\n',
'int m = rows(A);\n',
'vector[m] root_root_evals = sqrt(sqrt(eigenvalues_sym(A)));\n',
'matrix[m, m] evecs = eigenvectors_sym(A);\n',
'matrix[m, m] eprod = diag_post_multiply(evecs, root_root_evals);\n',
'return tcrossprod(eprod);\n',
'}\n',
'/* Function to transform P_real to P */\n',
'/* see Heaps 2022 for details (https://doi.org/10.1080/10618600.2022.2079648)*/\n',
'matrix P_realtoP(matrix P_real) {\n',
'int m = rows(P_real);\n',
'matrix[m, m] B = tcrossprod(P_real);\n',
'for(i in 1:m) B[i, i] += 1.0;\n',
'return mdivide_left_spd(sqrtm(B), P_real);\n',
'}\n',
'/* Function to perform the reverse mapping*/\n',
'/* see Heaps 2022 for details (https://doi.org/10.1080/10618600.2022.2079648)*/\n',
'matrix[,] rev_mapping(matrix[] P, matrix Sigma) {\n',
'int p = size(P);\n',
'int m = rows(Sigma);\n',
'matrix[m, m] phi_for[p, p]; matrix[m, m] phi_rev[p, p];\n',
'matrix[m, m] Sigma_for[p+1]; matrix[m, m] Sigma_rev[p+1];\n',
'matrix[m, m] S_for; matrix[m, m] S_rev;\n',
'matrix[m, m] S_for_list[p+1];\n',
'matrix[m, m] Gamma_trans[p+1];\n',
'matrix[m, m] phiGamma[2, p];\n',
'// Step 1:\n',
'Sigma_for[p+1] = Sigma;\n',
'S_for_list[p+1] = sqrtm(Sigma);\n',
'for(s in 1:p) {\n',
'// In this block of code S_rev is B^{-1} and S_for is a working matrix\n',
'S_for = - tcrossprod(P[p-s+1]);\n',
'for(i in 1:m) S_for[i, i] += 1.0;\n',
'S_rev = sqrtm(S_for);\n',
'S_for_list[p-s+1] = mdivide_right_spd(mdivide_left_spd(S_rev,\n',
'sqrtm(quad_form_sym(Sigma_for[p-s+2], S_rev))), S_rev);\n',
'Sigma_for[p-s+1] = tcrossprod(S_for_list[p-s+1]);\n',
'}\n',
'// Step 2:\n',
'Sigma_rev[1] = Sigma_for[1];\n',
'Gamma_trans[1] = Sigma_for[1];\n',
'for(s in 0:(p-1)) {\n',
'S_for = S_for_list[s+1];\n',
'S_rev = sqrtm(Sigma_rev[s+1]);\n',
'phi_for[s+1, s+1] = mdivide_right_spd(S_for * P[s+1], S_rev);\n',
"phi_rev[s+1, s+1] = mdivide_right_spd(S_rev * P[s+1]', S_for);\n",
'Gamma_trans[s+2] = phi_for[s+1, s+1] * Sigma_rev[s+1];\n',
'if(s>=1) {\n',
'for(k in 1:s) {\n',
'phi_for[s+1, k] = phi_for[s, k] - phi_for[s+1, s+1] * phi_rev[s, s-k+1];\n',
'phi_rev[s+1, k] = phi_rev[s, k] - phi_rev[s+1, s+1] * phi_for[s, s-k+1];\n',
'}\n',
'for(k in 1:s) Gamma_trans[s+2] = Gamma_trans[s+2] + phi_for[s, k] * Gamma_trans[s+2-k];\n',
'}\n',
"Sigma_rev[s+2] = Sigma_rev[s+1] - quad_form_sym(Sigma_for[s+1],phi_rev[s+1, s+1]');\n",
'}\n',
'for(i in 1:p) phiGamma[1, i] = phi_for[p, i];\n',
"for(i in 1:p) phiGamma[2, i] = Gamma_trans[i]';\n",
'return phiGamma;\n',
'}\n',
'}'
)
}
model_file <- readLines(textConnection(model_file), n = -1)
# Add transformed data lines for Heaps stationarity constraints
if (any(grepl('transformed data {', model_file, fixed = TRUE))) {
model_file[grep('transformed data {', model_file, fixed = TRUE)] <- paste0(
'transformed data {\n',
'vector[n_series] trend_zeros = rep_vector(0.0, n_series);\n',
'// exchangeable partial autocorrelation hyperparameters\n',
'// see Heaps 2022 for details (https://doi.org/10.1080/10618600.2022.2079648)\n',
'vector[2] es;\n',
'vector<lower=0>[2] fs;\n',
'vector<lower=0>[2] gs;\n',
'vector<lower=0>[2] hs;\n',
'es[1] = 0;\nes[2] = 0;\n',
'fs[1] = sqrt(0.455);\nfs[2] = sqrt(0.455);\n',
'gs[1] = 1.365;\ngs[2] = 1.365;\n',
'hs[1] = 0.071175;\nhs[2] = 0.071175;\n'
)
} else {
params_line <- min(which(grepl('parameters {', model_file, fixed = TRUE)))
model_file[params_line] <- paste0(
'transformed data {\n',
'vector[n_series] trend_zeros = rep_vector(0.0, n_series);\n',
'// exchangeable partial autocorrelation hyperparameters\n',
'// see Heaps 2022 for details (https://doi.org/10.1080/10618600.2022.2079648)\n',
'vector[2] es;\n',
'vector<lower=0>[2] fs;\n',
'vector<lower=0>[2] gs;\n',
'vector<lower=0>[2] hs;\n',
'es[1] = 0;\nes[2] = 0;\n',
'fs[1] = sqrt(0.455);\nfs[2] = sqrt(0.455);\n',
'gs[1] = 1.365;\ngs[2] = 1.365;\n',
'hs[1] = 0.071175;\nhs[2] = 0.071175;\n',
'}\n\nparameters {'
)
model_file <- readLines(textConnection(model_file), n = -1)
}
# Add parameters for real-valued partial autocorrelations
model_file[grep(
"vector<lower=0>[n_series] sigma;",
model_file,
fixed = TRUE
)] <- paste0(
"cholesky_factor_corr[n_series] L_Omega;\n",
"vector<lower=0>[n_series] sigma;\n\n",
'// unconstrained VAR1 partial autocorrelations\n',
'matrix[n_series, n_series] P_real;\n',
'// partial autocorrelation hyperparameters\n',
'vector[2] Pmu;\n',
'vector<lower=0>[2] Pomega;\n'
)
model_file <- readLines(textConnection(model_file), n = -1)
# Add transformed parameters for partial autocorrelation reverse mapping
model_file[grep(
'transformed parameters {',
model_file,
fixed = TRUE
)] <- paste0(
'transformed parameters {\n',
"// latent trend VAR1 autoregressive terms\n",
"matrix[n_series, n_series] A;\n",
'// LKJ form of covariance matrix\n',
"matrix[n_series, n_series] L_Sigma;\n",
'// computed error covariance matrix\n',
'cov_matrix[n_series] Sigma;\n',
'// initial trend covariance\n',
'cov_matrix[n_series] Gamma;'
)
model_file[
grep(
'trend[i, 1:n_series] = to_row_vector(trend_raw[i]);',
model_file,
fixed = TRUE
) +
1
] <- paste0(
'}\n\n',
'L_Sigma = diag_pre_multiply(sigma, L_Omega);\n',
'Sigma = multiply_lower_tri_self_transpose(L_Sigma);\n',
'// stationary VAR reparameterisation\n',
'{\n',
'matrix[n_series, n_series] P[1];\n',
'matrix[n_series, n_series] phiGamma[2, 1];\n',
'P[1] = P_realtoP(P_real);\n',
'phiGamma = rev_mapping(P, Sigma);\n',
'A = phiGamma[1, 1];\n',
'Gamma = phiGamma[2, 1];\n',
'}'
)
model_file <- readLines(textConnection(model_file), n = -1)
# Add priors for partial autocorrelations
model_file[grep("// trend means", model_file, fixed = TRUE)] <- paste0(
'// LKJ error correlation prior\n',
'L_Omega ~ lkj_corr_cholesky(2);\n',
'// partial autocorrelation hyperpriors\n',
'Pmu ~ normal(es, fs);\n',
'Pomega ~ gamma(gs, hs);\n',
'// unconstrained partial autocorrelations\n',
'diagonal(P_real) ~ normal(Pmu[1], 1 / sqrt(Pomega[1]));\n',
'for(i in 1:n_series) {\n',
'for(j in 1:n_series) {\n',
'if(i != j) P_real[i, j] ~ normal(Pmu[2], 1 / sqrt(Pomega[2]));\n',
'}\n',
'}\n',
"// trend means"
)
model_file <- readLines(textConnection(model_file), n = -1)
# Update prior for error SD parameters sigma
model_file[
grep(
"// priors for latent trend variance parameters",
model_file,
fixed = TRUE
) +
1
] <- paste0(
'sigma ~ inv_gamma(2.3693353, 0.7311319);'
)
model_file <- readLines(textConnection(model_file), n = -1)
# Update trend model to use multinormal
model_file[grep(
"trend_raw[1] ~ normal(0, sigma);",
model_file,
fixed = TRUE
)] <- "trend_raw[1] ~ multi_normal(trend_zeros, Gamma);"
model_file[grep(
"trend_raw[i] ~ normal(mu[i - 1], sigma);",
model_file,
fixed = TRUE
)] <- "trend_raw[i] ~ multi_normal_cholesky(mu[i - 1], L_Sigma);"
model_file <- gsub(
'contemporaneously uncorrelated',
'contemporaneously correlated',
model_file
)
model_file <- readLines(textConnection(model_file), n = -1)
# Return the updated model file
return(model_file)
}
R Package Documentation
Browse R Packages
We want your feedback!
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
Embedding an R snippet on your website
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.