Nothing
#' @import stats
#'
#' @title Compute log-likelihood of a GMVAR, StMVAR, and G-StMVAR models
#'
#' @description \code{loglikelihood_int} computes log-likelihoodof a GMVAR, StMVAR, and G-StMVAR models.
#'
#' @inheritParams in_paramspace_int
#' @param data a matrix or class \code{'ts'} object with \code{d>1} columns. Each column is taken to represent
#' a univariate time series. \code{NA} values are not supported.
#' @param p a positive integer specifying the autoregressive order of the model.
#' @param M \describe{
#' \item{For \strong{GMVAR} and \strong{StMVAR} models:}{a positive integer specifying the number of mixture components.}
#' \item{For \strong{G-StMVAR} models:}{a size (2x1) integer vector specifying the number of \emph{GMVAR type} components \code{M1}
#' in the first element and \emph{StMVAR type} components \code{M2} in the second element. The total number of mixture components
#' is \code{M=M1+M2}.}
#' }
#' @param params a real valued vector specifying the parameter values.
#' \describe{
#' \item{\strong{For unconstrained models:}}{
#' Should be size \eqn{((M(pd^2+d+d(d+1)/2+2)-M1-1)x1)} and have the form
#' \strong{\eqn{\theta}}\eqn{ = }(\strong{\eqn{\upsilon}}\eqn{_{1}},
#' ...,\strong{\eqn{\upsilon}}\eqn{_{M}}, \eqn{\alpha_{1},...,\alpha_{M-1},}\strong{\eqn{\nu}}\eqn{)}, where
#' \itemize{
#' \item \strong{\eqn{\upsilon}}\eqn{_{m}} \eqn{ = (\phi_{m,0},}\strong{\eqn{\phi}}\eqn{_{m}}\eqn{,\sigma_{m})}
#' \item \strong{\eqn{\phi}}\eqn{_{m}}\eqn{ = (vec(A_{m,1}),...,vec(A_{m,p})}
#' \item and \eqn{\sigma_{m} = vech(\Omega_{m})}, m=1,...,M,
#' \item \strong{\eqn{\nu}}\eqn{=(\nu_{M1+1},...,\nu_{M})}
#' \item \eqn{M1} is the number of GMVAR type regimes.
#' }
#' }
#' \item{\strong{For constrained models:}}{
#' Should be size \eqn{((M(d+d(d+1)/2+2)+q-M1-1)x1)} and have the form
#' \strong{\eqn{\theta}}\eqn{ = (\phi_{1,0},...,\phi_{M,0},}\strong{\eqn{\psi},}
#' \eqn{\sigma_{1},...,\sigma_{M},\alpha_{1},...,\alpha_{M-1},}\strong{\eqn{\nu}}), where
#' \itemize{
#' \item \strong{\eqn{\psi}} \eqn{(qx1)} satisfies (\strong{\eqn{\phi}}\eqn{_{1}}\eqn{,...,}
#' \strong{\eqn{\phi}}\eqn{_{M}) =} \strong{\eqn{C \psi}} where \strong{\eqn{C}} is a \eqn{(Mpd^2xq)}
#' constraint matrix.
#' }
#' }
#' \item{\strong{For same_means models:}}{
#' Should have the form
#' \strong{\eqn{\theta}}\eqn{ = (}\strong{\eqn{\mu},}\strong{\eqn{\psi},}
#' \eqn{\sigma_{1},...,\sigma_{M},\alpha_{1},...,\alpha_{M-1},}\strong{\eqn{\nu}}\eqn{)}, where
#' \itemize{
#' \item \strong{\eqn{\mu}}\eqn{= (\mu_{1},...,\mu_{g})} where
#' \eqn{\mu_{i}} is the mean parameter for group \eqn{i} and
#' \eqn{g} is the number of groups.
#' \item If AR constraints are employed, \strong{\eqn{\psi}} is as for constrained
#' models, and if AR constraints are not employed, \strong{\eqn{\psi}}\eqn{ = }
#' (\strong{\eqn{\phi}}\eqn{_{1}}\eqn{,...,}\strong{\eqn{\phi}}\eqn{_{M})}.
#' }
#' }
#' \item{\strong{For models with weight_constraints:}}{Drop \eqn{\alpha_1,...,\alpha_{M-1}} from
#' the parameter vector.}
#' \item{\strong{For structural models:}}{
#' Reduced form models can be directly used as recursively identified structural models. If the structural model is
#' identified by conditional heteroskedasticity, the parameter vector should have the form
#' \strong{\eqn{\theta}}\eqn{ = (\phi_{1,0},...,\phi_{M,0},}\strong{\eqn{\phi}}\eqn{_{1},...,}\strong{\eqn{\phi}}\eqn{_{M},
#' vec(W),}\strong{\eqn{\lambda}}\eqn{_{2},...,}\strong{\eqn{\lambda}}\eqn{_{M},\alpha_{1},...,\alpha_{M-1},}\strong{\eqn{\nu}}\eqn{)},
#' where
#' \itemize{
#' \item\strong{\eqn{\lambda}}\eqn{_{m}=(\lambda_{m1},...,\lambda_{md})} contains the eigenvalues of the \eqn{m}th mixture component.
#' }
#' \describe{
#' \item{\strong{If AR parameters are constrained: }}{Replace \strong{\eqn{\phi}}\eqn{_{1}}\eqn{,...,}
#' \strong{\eqn{\phi}}\eqn{_{M}} with \strong{\eqn{\psi}} \eqn{(qx1)} that satisfies (\strong{\eqn{\phi}}\eqn{_{1}}\eqn{,...,}
#' \strong{\eqn{\phi}}\eqn{_{M}) =} \strong{\eqn{C \psi}}, as above.}
#' \item{\strong{If same_means: }}{Replace \eqn{(\phi_{1,0},...,\phi_{M,0})} with \eqn{(\mu_{1},...,\mu_{g})},
#' as above.}
#' \item{\strong{If \eqn{W} is constrained:}}{Remove the zeros from \eqn{vec(W)} and make sure the other entries satisfy
#' the sign constraints.}
#' \item{\strong{If \eqn{\lambda_{mi}} are constrained via \code{C_lambda}:}}{Replace \strong{\eqn{\lambda}}\eqn{_{2},...,}
#' \strong{\eqn{\lambda}}\eqn{_{M}} with \strong{\eqn{\gamma}} \eqn{(rx1)} that satisfies (\strong{\eqn{\lambda}}\eqn{_{2}}
#' \eqn{,...,} \strong{\eqn{\lambda}}\eqn{_{M}) =} \strong{\eqn{C_{\lambda} \gamma}} where \eqn{C_{\lambda}} is
#' a \eqn{(d(M-1) x r)} constraint matrix.}
#' \item{\strong{If \eqn{\lambda_{mi}} are constrained via \code{fixed_lambdas}:}}{Drop \strong{\eqn{\lambda}}\eqn{_{2},...,}
#' \strong{\eqn{\lambda}}\eqn{_{M}} from the parameter vector.}
#' }
#' }
#' }
#' Above, \eqn{\phi_{m,0}} is the intercept parameter, \eqn{A_{m,i}} denotes the \eqn{i}th coefficient matrix of the \eqn{m}th
#' mixture component, \eqn{\Omega_{m}} denotes the error term covariance matrix of the \eqn{m}:th mixture component, and
#' \eqn{\alpha_{m}} is the mixing weight parameter. The \eqn{W} and \eqn{\lambda_{mi}} are structural parameters replacing the
#' error term covariance matrices (see Virolainen, 2022). If \eqn{M=1}, \eqn{\alpha_{m}} and \eqn{\lambda_{mi}} are dropped.
#' If \code{parametrization=="mean"}, just replace each \eqn{\phi_{m,0}} with regimewise mean \eqn{\mu_{m}}.
#' \eqn{vec()} is vectorization operator that stacks columns of a given matrix into a vector. \eqn{vech()} stacks columns
#' of a given matrix from the principal diagonal downwards (including elements on the diagonal) into a vector.
#'
#' In the \strong{GMVAR model}, \eqn{M1=M} and \strong{\eqn{\nu}} is dropped from the parameter vector. In the \strong{StMVAR} model,
#' \eqn{M1=0}. In the \strong{G-StMVAR} model, the first \code{M1} regimes are \emph{GMVAR type} and the rest \code{M2} regimes are
#' \emph{StMVAR type}. In \strong{StMVAR} and \strong{G-StMVAR} models, the degrees of freedom parameters in \strong{\eqn{\nu}} should
#' be strictly larger than two.
#'
#' The notation is similar to the cited literature.
#' @param model is "GMVAR", "StMVAR", or "G-StMVAR" model considered? In the G-StMVAR model, the first \code{M1} components
#' are GMVAR type and the rest \code{M2} components are StMVAR type.
#' @param conditional a logical argument specifying whether the conditional or exact log-likelihood function
#' should be used.
#' @param parametrization \code{"intercept"} or \code{"mean"} determining whether the model is parametrized with intercept
#' parameters \eqn{\phi_{m,0}} or regime means \eqn{\mu_{m}}, m=1,...,M.
#' @param constraints a size \eqn{(Mpd^2 x q)} constraint matrix \strong{\eqn{C}} specifying general linear constraints
#' to the autoregressive parameters. We consider constraints of form
#' (\strong{\eqn{\phi}}\eqn{_{1}}\eqn{,...,}\strong{\eqn{\phi}}\eqn{_{M}) = }\strong{\eqn{C \psi}},
#' where \strong{\eqn{\phi}}\eqn{_{m}}\eqn{ = (vec(A_{m,1}),...,vec(A_{m,p}) (pd^2 x 1), m=1,...,M},
#' contains the coefficient matrices and \strong{\eqn{\psi}} \eqn{(q x 1)} contains the related parameters.
#' For example, to restrict the AR-parameters to be the same for all regimes, set \strong{\eqn{C}}=
#' [\code{I:...:I}]\strong{'} \eqn{(Mpd^2 x pd^2)} where \code{I = diag(p*d^2)}.
#' Ignore (or set to \code{NULL}) if linear constraints should \strong{not} be employed.
#' @param same_means Restrict the mean parameters of some regimes to be the same? Provide a list of numeric vectors
#' such that each numeric vector contains the regimes that should share the common mean parameters. For instance, if
#' \code{M=3}, the argument \code{list(1, 2:3)} restricts the mean parameters of the second and third regime to be
#' the same but the first regime has freely estimated (unconditional) mean. Ignore or set to \code{NULL} if mean parameters
#' should not be restricted to be the same among any regimes. \strong{This constraint is available only for mean parametrized models;
#' that is, when \code{parametrization="mean"}.}
#' @param weight_constraints a numeric vector of length \eqn{M-1} specifying fixed parameter values for the mixing weight parameters
#' \eqn{\alpha_m, \ m=1,...,M-1}. Each element should be strictly between zero and one, and the sum of all the elements should
#' be strictly less than one.
#' @param structural_pars If \code{NULL} a reduced form model is considered. Reduced models can be used directly as recursively
#' identified structural models. For a structural model identified by conditional heteroskedasticity, should be a list containing
#' at least the first one of the following elements:
#' \itemize{
#' \item \code{W} - a \eqn{(dxd)} matrix with its entries imposing constraints on \eqn{W}: \code{NA} indicating that the element is
#' unconstrained, a positive value indicating strict positive sign constraint, a negative value indicating strict
#' negative sign constraint, and zero indicating that the element is constrained to zero.
#' \item \code{C_lambda} - a \eqn{(d(M-1) x r)} constraint matrix that satisfies (\strong{\eqn{\lambda}}\eqn{_{2}}\eqn{,...,}
#' \strong{\eqn{\lambda}}\eqn{_{M}) =} \strong{\eqn{C_{\lambda} \gamma}} where \strong{\eqn{\gamma}} is the new \eqn{(r x 1)}
#' parameter subject to which the model is estimated (similarly to AR parameter constraints). The entries of \code{C_lambda}
#' must be either \strong{positive} or \strong{zero}. Ignore (or set to \code{NULL}) if the eigenvalues \eqn{\lambda_{mi}}
#' should not be constrained.
#' \item \code{fixed_lambdas} - a length \eqn{d(M-1)} numeric vector (\strong{\eqn{\lambda}}\eqn{_{2}}\eqn{,...,}
#' \strong{\eqn{\lambda}}\eqn{_{M})} with elements strictly larger than zero specifying the fixed parameter values for the
#' parameters \eqn{\lambda_{mi}} should be constrained to. This constraint is alternative \code{C_lambda}.
#' Ignore (or set to \code{NULL}) if the eigenvalues \eqn{\lambda_{mi}} should not be constrained.
#' }
#' See Virolainen (forthcoming) for the conditions required to identify the shocks and for the B-matrix as well (it is \eqn{W} times
#' a time-varying diagonal matrix with positive diagonal entries).
#' @param check_params should it be checked that the parameter vector satisfies the model assumptions? Can be skipped to save
#' computation time if it does for sure.
#' @param minval the value that will be returned if the parameter vector does not lie in the parameter space
#' (excluding the identification condition).
#' @param to_return should the returned object be the log-likelihood value, which is default, or something else?
#' See the section "Return" for all the options.
#' @details \code{loglikelihood_int} takes use of the function \code{dmvn} from the package \code{mvnfast}.
#' @return
#' \describe{
#' \item{By default:}{log-likelihood value of the specified GMVAR, StMVAR, or G-StMVAR model,}
#' \item{If \code{to_return=="mw"}:}{a size ((n_obs-p)xM) matrix containing the mixing weights: for m:th component in m:th column.}
#' \item{If \code{to_return=="mw_tplus1"}:}{a size ((n_obs-p+1)xM) matrix containing the mixing weights: for m:th component in m:th column.
#' The last row is for \eqn{\alpha_{m,T+1}}.}
#' \item{If \code{to_return=="terms"}:}{a size ((n_obs-p)x1) numeric vector containing the terms \eqn{l_{t}}.}
#' \item{if \code{to_return=="loglik_and_mw"}:}{a list of two elements. The first element contains the log-likelihood value and the
#' second element contains the mixing weights.}
#' \item{If \code{to_return=="regime_cmeans"}:}{an \code{[T-p, d, M]} array containing the regimewise conditional means
#' (the first p values are used as the initial values).}
#' \item{If \code{to_return=="regime_ccovs"}:}{an \code{[d, d, T-p, M]} array containing the regimewise conditional
#' covariance matrices (the first p values are used as the initial values). The index \code{[ , , t, m]} gives the time
#' \code{t} conditional covariance matrix for the regime \code{m}.}
#' \item{If \code{to_return=="total_cmeans"}:}{a \code{[T-p, d]} matrix containing the conditional means of the process
#' (the first p values are used as the initial values).}
#' \item{If \code{to_return=="total_ccov"}:}{an \code{[d, d, T-p]} array containing the conditional covariance matrices of the process
#' (the first p values are used as the initial values).}
#' \item{If \code{to_return=="arch_scalars"}:}{a \code{[T-p, M]} matrix containing the regimewise arch scalars
#' multiplying error term covariance matrix in the conditional covariance matrix of the regime. For GMVAR type regimes, these
#' are all ones (the first p values are used as the initial values).}
#' \item{if \code{to_return=="loglik_mw_archscalars"}:}{a list of three elements. The first element contains the log-likelihood value, the
#' second element contains the mixing weights, the third element contains the arch scalars
#' (this is used in \code{quantile_residuals_int}).}
#' }
#' @references
#' \itemize{
#' \item Kalliovirta L., Meitz M. and Saikkonen P. 2016. Gaussian mixture vector autoregression.
#' \emph{Journal of Econometrics}, \strong{192}, 485-498.
#' \item Lütkepohl H. 2005. New Introduction to Multiple Time Series Analysis,
#' \emph{Springer}.
#' \item McElroy T. 2017. Computation of vector ARMA autocovariances.
#' \emph{Statistics and Probability Letters}, \strong{124}, 92-96.
#' \item Virolainen S. (forthcoming). A statistically identified structural vector autoregression with endogenously
#' switching volatility regime. \emph{Journal of Business & Economic Statistics}.
#' \item Virolainen S. 2022. Gaussian and Student's t mixture vector autoregressive model with application to the
#' asymmetric effects of monetary policy shocks in the Euro area. Unpublished working
#' paper, available as arXiv:2109.13648.
#' }
#' @keywords internal
loglikelihood_int <- function(data, p, M, params, model=c("GMVAR", "StMVAR", "G-StMVAR"), conditional=TRUE,
parametrization=c("intercept", "mean"), constraints=NULL, same_means=NULL,
weight_constraints=NULL, structural_pars=NULL,
to_return=c("loglik", "mw", "mw_tplus1", "loglik_and_mw", "terms",
"regime_cmeans", "regime_ccovs", "total_cmeans", "total_ccovs",
"arch_scalars", "loglik_mw_archscalars"),
check_params=TRUE, minval=NULL, stat_tol=1e-3, posdef_tol=1e-8, df_tol=1e-8) {
# Compute required values
epsilon <- round(log(.Machine$double.xmin) + 10) # Logarithm of the smallest value that can be handled normally
d <- ncol(data)
n_obs <- nrow(data)
T_obs <- n_obs - p
to_return <- match.arg(to_return)
# Collect parameter values
model <- match.arg(model)
parametrization <- match.arg(parametrization)
check_same_means(parametrization=parametrization, same_means=same_means)
params <- reform_constrained_pars(p=p, M=M, d=d, params=params, model=model, constraints=constraints,
same_means=same_means, weight_constraints=weight_constraints,
structural_pars=structural_pars) # All constraints are expanded and removed from the parameter vector
W_constraints <- structural_pars$W
structural_pars <- get_unconstrained_structural_pars(structural_pars=structural_pars)
if(parametrization == "intercept") {
all_phi0 <- pick_phi0(p=p, M=M, d=d, params=params, structural_pars=structural_pars)
} else {
mu <- pick_phi0(p=p, M=M, d=d, params=params, structural_pars=structural_pars) # mean parameters instead of phi0
}
all_A <- pick_allA(p=p, M=M, d=d, params=params, structural_pars=structural_pars) # A_{m,i}, m=1,...,M, i=1,..,p
all_Omega <- pick_Omegas(p=p, M=M, d=d, params=params, structural_pars=structural_pars) # Omega_m
all_boldA <- form_boldA(p=p, M=M, d=d, all_A=all_A) # The 'bold A' for each m=1,..,M, Lütkepohl 2005, eq.(2.1.8)
alphas <- pick_alphas(p=p, M=M, d=d, params=params, model=model) # Mixing weight parameters
all_df <- pick_df(M=M, params=params, model=model) # Degrees of freedom parameters (numeric(0) for GMVAR models)
# Check that the parameter vector lies in the parameter space (excluding indentifiability)
if(check_params) {
if(!in_paramspace_int(p=p, M=M, d=d, params=params, model=model, all_boldA=all_boldA, alphas=alphas, all_Omega=all_Omega,
W_constraints=W_constraints, stat_tol=stat_tol, posdef_tol=posdef_tol, df_tol=df_tol)) {
return(minval)
}
}
if(model == "GMVAR") {
M1 <- M
M2 <- 0
} else if(model == "StMVAR") {
M1 <- 0
M2 <- M
} else { # model == "G-StMVAR"
M1 <- M[1]
M2 <- M[2]
}
M <- sum(M) # The total number of mixture components
# i:th row denotes the vector \bold{y_{i-1}} = (y_{i-1},...,y_{i-p}) (dpx1),
# assuming the observed data is y_{-p+1},...,y_0,y_1,...,y_{T}
Y <- reform_data(data, p)
# Calculate expected values (column per component) or phi0-parameters if using mean-parametrization
Id <- diag(nrow=d)
if(parametrization == "intercept") {
mu <- vapply(1:M, function(m) solve(Id - rowSums(all_A[, , , m, drop=FALSE], dims=2), all_phi0[,m]), numeric(d)) # rowSums: sum over dims+1=3
} else {
all_phi0 <- vapply(1:M, function(m) (Id - rowSums(all_A[, , , m, drop=FALSE], dims=2))%*%mu[,m], numeric(d))
}
# Calculate the covariance matrices Sigma_{m,p} (Lutkepohl 2005, eq. (2.1.39) or the algorithm proposed by McElroy 2017)
Sigmas <- get_Sigmas(p=p, M=M, d=d, all_A=all_A, all_boldA=all_boldA, all_Omega=all_Omega) # Store the (dpxdp) covariance matrices
inv_Sigmas <- array(dim=c(d*p, d*p, M)) # Only used for StMVAR type regimes
chol_Sigmas <- array(dim=c(d*p, d*p, M))
for(m in 1:M) {
# Take Cholesky here to avoid unnecessary warnings from mvnfast::dmvn, also used in inverting for m > M1:
chol_Sigmas[, , m] <- chol(Sigmas[, , m])
if(m > M1) {
inv_Sigmas[, , m] <- chol2inv(chol_Sigmas[, , m])
}
}
# Calculate the dp-dimensional multinormal densities (KMS 2016, eq.(6)) or log Students t densities (Virolainen 2022, eq. (3.4)):
# i:th row for index i-1 etc, m:th column for m:th component.
# We calculate in logarithm because the non-log values may be too close to zero for machine accuracy (if they are too close to zero
# for all regimes and computer handles them as zero, we would divide by zero when calculating the mixing weights)
log_mvdvalues <- matprods <- matrix(nrow=n_obs - p + 1, ncol=M) # The quadratic forms in Student's t density
if(M1 > 0) { # Multinormals
log_mvdvalues[,1:M1] <- vapply(1:M1, function(m) mvnfast::dmvn(X=Y, mu=rep(mu[,m], p), sigma=chol_Sigmas[, , m],
log=TRUE, ncores=1, isChol=TRUE), numeric(T_obs + 1))
}
if(M2 > 0) { # Multistudents
for(m in (M1 + 1):M) { # Go through the StMVAR type regimes
tmp_mat <- Y - matrix(rep(mu[,m], p), nrow=nrow(Y), ncol=ncol(Y), byrow=TRUE)
matprods[,m] <- rowSums(tmp_mat%*%inv_Sigmas[, , m]*tmp_mat)
logC <- lgamma(0.5*(d*p + all_df[m - M1])) - 0.5*d*p*log(base::pi) - 0.5*d*p*log(all_df[m - M1] - 2) - lgamma(0.5*all_df[m - M1])
log_det_Sigma <- 2*log(prod(diag(chol_Sigmas[, , m]))) #log(det(Sigmas[, , m]))
log_mvdvalues[,m] <- logC - 0.5*log_det_Sigma - 0.5*(d*p + all_df[m - M1])*log(1 + matprods[,m]/(all_df[m - M1] - 2))
}
}
## Calculate the mixing weights alpha_{m,t} (KMS 2016, eq.(7))
if(to_return != "mw_tplus1") {
log_mvdvalues <- log_mvdvalues[1:T_obs, , drop=FALSE] # alpha_mt uses y_{t-1} so the last row is not needed
}
alpha_mt_and_l_0 <- get_alpha_mt(M=M, log_mvdvalues=log_mvdvalues, alphas=alphas,
epsilon=epsilon, conditional=conditional, also_l_0=TRUE)
alpha_mt <- alpha_mt_and_l_0$alpha_mt
l_0 <- alpha_mt_and_l_0$l_0 # The first term in the exact log-likelihood function (=0 for conditional)
if(to_return == "mw" | to_return == "mw_tplus1") {
return(alpha_mt)
}
# Calculate the conditional means mu_{m,t} (KMS 2016, Condition 1 (a))
# The dimensions of mu_mt will be: [t, p, m]
all_A2 <- array(all_A, dim=c(d, d*p, M)) # cbind coefficient matrices of each component: m:th component is obtained at [, , m]
Y2 <- Y[1:T_obs,] # Last row is not needed because mu_mt uses lagged values
mu_mt <- array(vapply(1:M, function(m) t(all_phi0[, m] + tcrossprod(all_A2[, , m], Y2)), numeric(d*T_obs)), dim=c(T_obs, d, M)) # [, , m]
# For the StMVAR type regimes, calculate the time-varying ARC-type scalar that multiplies the conditional covariance matrix:
# i:th row calculated from y_{i-1} observation when the indexing is y_{-p+1},..,y_0,y_1,...,y_T - note that the last row of
# arch_scalars is for time T+1
arch_scalars <- matrix(1, nrow=nrow(matprods) - 1, ncol=M) # The arch-scalars are all 1 for the GMVAR type regimes
if(M2 > 0) {
# The last row is not needed because for the time t, we use the previous value y_{t-1}:
matprods0 <- as.matrix(matprods[1:(nrow(matprods) - 1), (M1 + 1):M])
tmp_mat <- matrix(all_df, nrow=nrow(matprods0), ncol=ncol(matprods0), byrow=TRUE)
arch_scalars[,(M1 + 1):M] <- (tmp_mat - 2 + matprods0)/(tmp_mat - 2 + d*p)
}
if(to_return == "regime_cmeans") {
return(mu_mt)
} else if(to_return == "arch_scalars") {
return(arch_scalars)
} else if(to_return == "total_cmeans") { # KMS 2016, eq.(3)
return(matrix(rowSums(vapply(1:M, function(m) alpha_mt[,m]*mu_mt[, , m], numeric(d*T_obs))), nrow=T_obs, ncol=d, byrow=FALSE))
} else if(to_return == "regime_ccovs") {
regime_ccovs <- array(vapply(1:M, function(m) t(matrix(all_Omega[, , m], nrow=nrow(arch_scalars), ncol=d^2, byrow=TRUE)*arch_scalars[, m]),
numeric(d*d*nrow(arch_scalars))), dim=c(d, d, nrow(arch_scalars), M))
return(regime_ccovs)
} else if(to_return == "total_ccovs") { # KMS 2016, eq.(4)
first_term <- array(rowSums(vapply(1:M, function(m) rep(alpha_mt[, m]*arch_scalars[, m], each=d*d)*as.vector(all_Omega[, , m]),
numeric(d*d*T_obs))), dim=c(d, d, T_obs))
sum_alpha_mu <- matrix(rowSums(vapply(1:M, function(m) alpha_mt[, m]*mu_mt[, , m],
numeric(d*T_obs))), nrow=T_obs, ncol=d, byrow=FALSE)
second_term <- array(rowSums(vapply(1:M, function(m) rep(alpha_mt[, m],
each=d*d)*as.vector(vapply(1:nrow(alpha_mt),
function(i1) tcrossprod((mu_mt[, , m] - sum_alpha_mu)[i1,]),
numeric(d*d))), numeric(d*d*T_obs))), dim=c(d, d, T_obs))
return(first_term + second_term)
}
# Calculate the second term of the log-likelihood (KMS 2016 eq.(10)), also see Virolainen (2022) for the StMVAR type regimes
dat <- data[(p + 1):n_obs,] # Initial values are not used here (conditional means and variances are already calculated)
mvd_vals <- matrix(nrow=T_obs, ncol=M)
if(M1 > 0) { # GMVAR type regimes
mvd_vals[,1:M1] <- vapply(1:M1, function(m) mvnfast::dmvn(X=dat - mu_mt[, , m], mu=rep(0, times=d), sigma=all_Omega[, , m],
log=FALSE, ncores=1, isChol=FALSE), numeric(T_obs))
}
if(M2 > 0) { # StMVAR type regimes
for(m in (M1 + 1):M) {
df_m <- all_df[m - M1] + d*p # Degrees of freedom in regime m
chol_Omega_m <- chol(all_Omega[, , m]) # Faster determinant and matrix inversion useing Cholesky decomposition
det_Omega_m <- prod(diag(chol_Omega_m))^2
inv_Omega_m <- chol2inv(chol_Omega_m)
# Below, we calculate the d-dimensional conditional t-densities for the regime m, for t=1,...,T.
tmp_mat <- dat - mu_mt[, , m]
mvd_vals[, m] <- exp(lgamma(0.5*(d + df_m)) - 0.5*d*log(base::pi) - 0.5*d*log(df_m - 2) - lgamma(0.5*df_m))*
arch_scalars[, m]^(-d/2)*1/sqrt(det_Omega_m)*(1 + rowSums(tmp_mat%*%inv_Omega_m*tmp_mat)/(arch_scalars[, m]*(df_m - 2)))^(-0.5*(d + df_m))
}
}
weighted_mvd <- rowSums(alpha_mt*mvd_vals)
weighted_mvd[weighted_mvd == 0] <- exp(epsilon)
l_t <- log(weighted_mvd)
loglik <- l_0 + sum(l_t)
if(to_return == "terms") {
return(l_t)
} else if(to_return == "loglik_and_mw") {
return(list(loglik=loglik,
mw=alpha_mt))
} else if(to_return == "loglik_mw_archscalars") {
return(list(loglik=loglik,
mw=alpha_mt,
arch_scalars=arch_scalars))
} else { # to_return == "loglik"
return(loglik)
}
}
#' @title Get mixing weights alpha_mt (this function is for internal use)
#'
#' @description \code{get_alpha_mt} computes the mixing weights based on
#' the logarithm of the multivariate normal densities in the definition of
#' the mixing weights.
#'
#' @inheritParams loglikelihood_int
#' @param log_mvdvalues \eqn{T x M} matrix containing the log multivariate normal densities.
#' @param alphas \eqn{M x 1} vector containing the mixing weight pa
#' @param epsilon the smallest number such that its exponent is wont classified as numerically zero
#' (around \code{-698} is used).
#' @param also_l_0 return also l_0 (the first term in the exact log-likelihood function)?
#' @details Note that we index the time series as \eqn{-p+1,...,0,1,...,T} as in Kalliovirta et al. (2016).
#' @return Returns the mixing weights a matrix of the same dimension as \code{log_mvdvalues} so
#' that the t:th row is for the time point t and m:th column is for the regime m.
#' @inherit in_paramspace_int references
#' @seealso \code{\link{loglikelihood_int}}
#' @keywords internal
get_alpha_mt <- function(M, log_mvdvalues, alphas, epsilon, conditional, also_l_0=FALSE) {
if(M == 1) {
if(!is.matrix(log_mvdvalues)) log_mvdvalues <- as.matrix(log_mvdvalues) # Possibly many time points but only one regime
alpha_mt <- as.matrix(rep(1, nrow(log_mvdvalues)))
} else {
if(!is.matrix(log_mvdvalues)) log_mvdvalues <- t(as.matrix(log_mvdvalues)) # Only one time point but multiple regimes
log_mvdvalues_orig <- log_mvdvalues
small_logmvns <- log_mvdvalues < epsilon
if(any(small_logmvns)) {
# If too small or large non-log-density values are present (i.e., that would yield -Inf or Inf),
# we replace them with ones that are not too small or large but imply the same mixing weights
# up to negligible numerical tolerance.
which_change <- rowSums(small_logmvns) > 0 # Which rows contain too small values
to_change <- log_mvdvalues[which_change, , drop=FALSE]
largest_vals <- do.call(pmax, split(to_change, f=rep(1:ncol(to_change), each=nrow(to_change)))) # The largest values of those rows
diff_to_largest <- to_change - largest_vals # Differences to the largest value of the row
# For each element in each row, check the (negative) distance from the largest value of the row. If the difference
# is smaller than epsilon, replace the with epsilon. The results are then the new log_mvn values.
diff_to_largest[diff_to_largest < epsilon] <- epsilon
# Replace the old log_mvdvalues with the new ones
log_mvdvalues[which_change,] <- diff_to_largest
}
mvnvalues <- exp(log_mvdvalues)
denominator <- as.vector(mvnvalues%*%alphas)
alpha_mt <- (mvnvalues/denominator)%*%diag(alphas)
}
if(!also_l_0) {
return(alpha_mt)
} else {
# First term of the exact log-likelihood (Kalliovirta et al. 2016, eq.(9))
l_0 <- 0
if(M == 1 && conditional == FALSE) {
l_0 <- log_mvdvalues[1,]
} else if(M > 1 && conditional == FALSE) {
if(any(log_mvdvalues_orig[1,] < epsilon)) { # Need to use Brobdingnag
l_0 <- log(Reduce("+", lapply(1:M, function(i1) alphas[i1]*exp(Brobdingnag::as.brob(log_mvdvalues_orig[1, i1])))))
} else {
l_0 <- log(sum(alphas*mvnvalues[1,]))
}
}
return(list(alpha_mt=alpha_mt,
l_0=l_0))
}
}
#' @title Compute log-likelihood of a GMVAR, StMVAR, or G-StMVAR model using parameter vector
#'
#' @description \code{loglikelihood} computes log-likelihood of a GMVAR, StMVAR, or G-StMVAR model using parameter vector
#' instead of an object of class 'gsmvar'. Exists for convenience if one wants to for example
#' employ other estimation algorithms than the ones used in \code{fitGSMVAR}. Use \code{minval} to
#' control what happens when the parameter vector is outside the parameter space.
#'
#' @inheritParams loglikelihood_int
#' @return Returns log-likelihood if \code{params} is in the parameters space and \code{minval} if not.
#' @inherit loglikelihood_int details references
#' @seealso \code{\link{fitGSMVAR}}, \code{\link{GSMVAR}}, \code{\link{calc_gradient}}
#' @examples
#' # GMVAR(2, 2), d=2 model;
#' params22 <- c(0.36, 0.121, 0.223, 0.059, -0.151, 0.395, 0.406, -0.005,
#' 0.083, 0.299, 0.215, 0.002, 0.03, 0.484, 0.072, 0.218, 0.02, -0.119,
#' 0.722, 0.093, 0.032, 0.044, 0.191, 1.101, -0.004, 0.105, 0.58)
#' loglikelihood(data=gdpdef, p=2, M=2, params=params22)
#'
#' # Structural GMVAR(2, 2), d=2 model identified with sign-constraints:
#' params22s <- c(0.36, 0.121, 0.484, 0.072, 0.223, 0.059, -0.151, 0.395,
#' 0.406, -0.005, 0.083, 0.299, 0.218, 0.02, -0.119, 0.722, 0.093, 0.032,
#' 0.044, 0.191, 0.057, 0.172, -0.46, 0.016, 3.518, 5.154, 0.58)
#' W_22 <- matrix(c(1, 1, -1, 1), nrow=2, byrow=FALSE)
#' loglikelihood(data=gdpdef, p=2, M=2, params=params22s, structural_pars=list(W=W_22))
#' @export
loglikelihood <- function(data, p, M, params, model=c("GMVAR", "StMVAR", "G-StMVAR"), conditional=TRUE, parametrization=c("intercept", "mean"),
constraints=NULL, same_means=NULL, weight_constraints=NULL, structural_pars=NULL, minval=NA,
stat_tol=1e-3, posdef_tol=1e-8, df_tol=1e-8) {
parametrization <- match.arg(parametrization)
model <- match.arg(model)
check_same_means(parametrization=parametrization, same_means=same_means)
data <- check_data(data, p)
d <- ncol(data)
check_pMd(p=p, M=M, d=d, model=model)
check_constraints(p=p, M=M, d=d, constraints=constraints, same_means=same_means, weight_constraints=weight_constraints,
structural_pars=structural_pars)
if(length(params) != n_params(p=p, M=M, d=d, model=model, constraints=constraints, same_means=same_means,
weight_constraints=weight_constraints, structural_pars=structural_pars)) {
stop("Parameter vector has wrong dimension")
}
loglikelihood_int(data=data, p=p, M=M, params=params, model=model, conditional=conditional, parametrization=parametrization,
constraints=constraints, same_means=same_means, weight_constraints=weight_constraints,
structural_pars=structural_pars, to_return="loglik", check_params=TRUE, minval=minval,
stat_tol=stat_tol, posdef_tol=posdef_tol, df_tol=df_tol)
}
#' @title Compute conditional moments of a GMVAR, StMVAR, or G-StMVAR model
#'
#' @description \code{cond_moments} compute conditional regimewise means, conditional means, and conditional covariance matrices
#' of a GMVAR, StMVAR, or G-StMVAR model.
#'
#' @inheritParams loglikelihood_int
#' @param to_return should the regimewise conditional means, total conditional means, or total conditional covariance matrices
#' be returned?
#' @details The first p values are used as the initial values, and by conditional we mean conditioning on the past. Formulas
#' for the conditional means and covariance matrices are given in equations (3) and (4) of KMS (2016).
#' @return
#' \describe{
#' \item{If \code{to_return=="regime_cmeans"}:}{an \code{[T-p, d, M]} array containing the regimewise conditional means
#' (the first p values are used as the initial values).}
#' \item{If \code{to_return=="regime_ccovs"}:}{an \code{[d, d, T-p, M]} array containing the regimewise conditional
#' covariance matrices (the first p values are used as the initial values). The index \code{[ , , t, m]} gives the time
#' \code{t} conditional covariance matrix for the regime \code{m}.}
#' \item{If \code{to_return=="total_cmeans"}:}{a \code{[T-p, d]} matrix containing the conditional means of the process
#' (the first p values are used as the initial values).}
#' \item{If \code{to_return=="total_ccov"}:}{an \code{[d, d, T-p]} array containing the conditional covariance matrices of the process
#' (the first p values are used as the initial values).}
#' \item{If \code{to_return=="arch_scalars"}:}{a \code{[T-p, M]} matrix containing the regimewise arch scalars
#' multiplying error term covariance matrix in the conditional covariance matrix of the regime. For GMVAR type regimes, these
#' are all ones (the first p values are used as the initial values).}
#' }
#' @inherit loglikelihood_int references
#' @family moment functions
#' @examples
#' # GMVAR(2, 2), d=2 model;
#' params22 <- c(0.36, 0.121, 0.223, 0.059, -0.151, 0.395, 0.406, -0.005,
#' 0.083, 0.299, 0.215, 0.002, 0.03, 0.484, 0.072, 0.218, 0.02, -0.119,
#' 0.722, 0.093, 0.032, 0.044, 0.191, 1.101, -0.004, 0.105, 0.58)
#' cond_moments(data=gdpdef, p=2, M=2, params=params22, to_return="regime_cmeans")
#' cond_moments(data=gdpdef, p=2, M=2, params=params22, to_return="total_cmeans")
#' cond_moments(data=gdpdef, p=2, M=2, params=params22, to_return="total_ccovs")
#' @export
cond_moments <- function(data, p, M, params, model=c("GMVAR", "StMVAR", "G-StMVAR"), parametrization=c("intercept", "mean"),
constraints=NULL, same_means=NULL, weight_constraints=NULL, structural_pars=NULL,
to_return=c("regime_cmeans", "regime_ccovs", "total_cmeans", "total_ccovs", "arch_scalars"),
minval=NA, stat_tol=1e-3, posdef_tol=1e-8, df_tol=1e-8) {
parametrization <- match.arg(parametrization)
model <- match.arg(model)
to_return <- match.arg(to_return)
check_same_means(parametrization=parametrization, same_means=same_means)
data <- check_data(data, p)
d <- ncol(data)
check_pMd(p=p, M=M, d=d, model=model)
check_constraints(p=p, M=M, d=d, constraints=constraints, same_means=same_means,
weight_constraints=weight_constraints, structural_pars=structural_pars)
if(length(params) != n_params(p=p, M=M, d=d, model=model, constraints=constraints, same_means=same_means,
weight_constraints=weight_constraints, structural_pars=structural_pars)) {
stop("Parameter vector has wrong dimension")
}
loglikelihood_int(data=data, p=p, M=M, params=params, model=model, parametrization=parametrization,
constraints=constraints, same_means=same_means, weight_constraints=weight_constraints,
structural_pars=structural_pars, to_return=to_return, check_params=TRUE, minval=minval,
stat_tol=stat_tol, posdef_tol=posdef_tol, df_tol=df_tol)
}
#' @title Calculate AIC, HQIC, and BIC
#'
#' @description \code{get_IC} calculates the information criteria values AIC, HQIC, and BIC.
#'
#' @param loglik log-likelihood value
#' @param npars number of (freely estimated) parameters in the model
#' @param obs numbers of observations with starting values excluded for conditional models.
#' @details Note that for conditional models with different autoregressive order p the
#' information criteria values are \strong{NOT} comparable.
#' @return Returns a data frame containing the information criteria values.
#' @keywords internal
get_IC <- function(loglik, npars, obs) {
AIC <- -2*loglik + 2*npars
HQIC <- -2*loglik + 2*npars*log(log(obs))
BIC <- -2*loglik + npars*log(obs)
data.frame(AIC=AIC, HQIC=HQIC, BIC=BIC)
}
#' @title Calculate multivariate Pearson residuals of a GMVAR, StMVAR, or G-StMVAR model
#'
#' @description \code{Pearson_residuals} calculates multivariate Pearson residuals for a GMVAR, StMVAR, or G-StMVAR model.
#'
#' @inheritParams quantile_residual_tests
#' @param standardize Should the residuals be standardized? Use \code{FALSE} to obtain raw residuals.
#' @return Returns \eqn{((n_obs-p) x d)} matrix containing the residuals,
#' \eqn{j}:th column corresponds to the time series in the \eqn{j}:th column of the data.
#' @inherit GSMVAR references
#' @seealso \code{\link{fitGSMVAR}}, \code{\link{GSMVAR}}, \code{\link{quantile_residuals}},
#' \code{\link{diagnostic_plot}}
#' @examples
#' # GMVAR(1,2), d=2 model:
#' params12 <- c(0.55, 0.112, 0.344, 0.055, -0.009, 0.718, 0.319, 0.005, 0.03,
#' 0.619, 0.173, 0.255, 0.017, -0.136, 0.858, 1.185, -0.012, 0.136, 0.674)
#' mod12 <- GSMVAR(gdpdef, p=1, M=2, params=params12)
#' Pearson_residuals(mod12, standardize=FALSE) # Raw residuals
#' Pearson_residuals(mod12, standardize=TRUE) # Standardized to identity cov.matrix.
#'
#' # Structural GMVAR(2, 2), d=2 model identified with sign-constraints:
#' params22s <- c(0.36, 0.121, 0.484, 0.072, 0.223, 0.059, -0.151, 0.395,
#' 0.406, -0.005, 0.083, 0.299, 0.218, 0.02, -0.119, 0.722, 0.093, 0.032,
#' 0.044, 0.191, 0.057, 0.172, -0.46, 0.016, 3.518, 5.154, 0.58)
#' W_22 <- matrix(c(1, 1, -1, 1), nrow=2, byrow=FALSE)
#' mod22s <- GSMVAR(gdpdef, p=2, M=2, params=params22s, structural_pars=list(W=W_22))
#' Pearson_residuals(mod22s, standardize=FALSE) # Raw residuals
#' Pearson_residuals(mod22s, standardize=TRUE) # Standardized to identity cov.matrix.
#' @export
Pearson_residuals <- function(gsmvar, standardize=TRUE) {
# Checks, etc
gsmvar <- gmvar_to_gsmvar(gsmvar) # Backward compatibility
check_gsmvar(gsmvar)
p <- gsmvar$model$p
M <- gsmvar$model$M
d <- gsmvar$model$d
data <- gsmvar$data
n_obs <- nrow(data)
T_obs <- n_obs - p
# Conditional means
mu_t <- loglikelihood_int(data=data, p=gsmvar$model$p, M=gsmvar$model$M, params=gsmvar$params,
model=gsmvar$model$model, conditional=gsmvar$model$conditional,
parametrization=gsmvar$model$parametrization,
constraints=gsmvar$model$constraints,
same_means=gsmvar$model$same_means,
weight_constraints=gsmvar$model$weight_constraints,
structural_pars=gsmvar$model$structural_pars,
to_return="total_cmeans", check_params=TRUE,
stat_tol=gsmvar$num_tols$stat_tol,
posdef_tol=gsmvar$num_tols$posdef_tol,
df_tol=gsmvar$num_tols$df_tol)
# Nonstandardized residuals
y_minus_mu <- data[(p + 1):nrow(data),] - mu_t # [T_obs, d]
if(!standardize) {
return(y_minus_mu)
}
# Conditional covariance matrices
Omega_t <- loglikelihood_int(data=data, p=gsmvar$model$p, M=gsmvar$model$M, params=gsmvar$params,
model=gsmvar$model$model, conditional=gsmvar$model$conditional,
parametrization=gsmvar$model$parametrization,
constraints=gsmvar$model$constraints,
same_means=gsmvar$model$same_means,
weight_constraints=gsmvar$model$weight_constraints,
structural_pars=gsmvar$model$structural_pars,
to_return="total_ccovs", check_params=TRUE,
stat_tol=gsmvar$num_tols$stat_tol,
posdef_tol=gsmvar$num_tols$posdef_tol,
df_tol=gsmvar$num_tols$df_tol)
all_residuals <- matrix(nrow=T_obs, ncol=d)
# Calculate the Pearson residuals
for(i1 in 1:T_obs) {
all_residuals[i1,] <- solve(unvec(d=d, a=get_symmetric_sqrt(Omega_t[, , i1])), y_minus_mu[i1,])
}
all_residuals
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.