R/hypoTests.R

In sstvars: Toolkit for Reduced Form and Structural Smooth Transition Vector Autoregressive Models

#' @title Perform Wald test for a STVAR model
#'
#' @description \code{Wald_test} performs a Wald test for a STVAR model
#'
#' @param stvar an object of class \code{'stvar'} generated by \code{fitSTVAR} or \code{STVAR}, containing
#'   the model specified by the alternative hypothesis (i.e., \strong{the unconstrained model}).
#' @param A a size \eqn{(k\times}\code{n_params}\eqn{)} matrix with full row rank specifying a part of the null hypothesis
#'   where \code{n_params} is the number of parameters in the (unconstrained) model.
#'   See details for more information.
#' @param c a length \eqn{k} vector specifying a part of the null hypothesis. See details for more information.
#' @details Denoting the true parameter value by \eqn{\theta_{0}}, we test the null hypothesis \eqn{A\theta_{0}=c}.
#'   Under the null, the test statistic is asymptotically \eqn{\chi^2}-distributed with \eqn{k}
#'   (\code{=nrow(A)}) degrees of freedom. The parameter \eqn{\theta_{0}} is assumed to have the same form as in
#'   the model supplied in the argument \code{stvar} and it is presented in the documentation of the argument
#'   \code{params} in the function \code{STVAR} (see \code{?STVAR}).
#'
#'   \strong{The test is based on the assumption of the standard result of asymptotic normality!}
#'   Also note that this function does \strong{not} check whether the model assumptions hold under the null.
#' @return A list with class "hypotest" containing the test results and arguments used to calculate the test.
#' @seealso \code{\link{LR_test}}, \code{\link{Rao_test}}, \code{\link{fitSTVAR}}, \code{\link{STVAR}},
#'   \code{\link{diagnostic_plot}}, \code{\link{profile_logliks}}, \code{\link{Portmanteau_test}}
#' @references
#'  \itemize{
#'    \item Buse A. (1982). The Likelihood Ratio, Wald, and Lagrange Multiplier Tests: An Expository Note.
#'      \emph{The American Statistician}, 36(3a), 153-157.
#'  }
#' @examples
#'  # Logistic Student's t STVAR with p=1, M=2, and the first lag of the second variable
#'  # as the switching variable (parameter values were obtained by maximum likelihood estimation;
#'  # fitSTVAR is not used here because the estimation is computationally demanding).
#'  params12 <- c(0.62906848, 0.14245295, 2.41245785, 0.66719269, 0.3534745, 0.06041779, -0.34909745,
#'   0.61783824, 0.125769, -0.04094521, -0.99122586, 0.63805416, 0.371575, 0.00314754, 0.03440824,
#'   1.29072533, -0.06067807, 0.18737385, 1.21813844, 5.00884263, 7.70111672)
#'  fit12 <- STVAR(data=gdpdef, p=1, M=2, params=params12, weight_function="logistic",
#'   weightfun_pars=c(2, 1), cond_dist="Student")
#'  fit12
#'
#'  # Test whether the location parameter equals 1.
#'  # For this model, the parameter vector has the length 21 and
#'  # location parameter is in the 19th element:
#'  A <- matrix(c(rep(0, times=18), 1, 0, 0), nrow=1, ncol=21)
#'  c <- 1
#'  Wald_test(fit12, A=A, c=c)
#'
#'  # Test whether the intercepts and autoregressive matrices are identical across the regimes:
#'  # fit12 has parameter vector of length 21. In the first regime, the intercepts are in the
#'  # elements 1,2 and the AR parameters in the elements 5,...,8. In the second regime,
#'  # the intercepts are in the elements 3,4, and the AR parameters the elements 9,...,12.
#'  A <- rbind(cbind(diag(2), -diag(2), matrix(0, nrow=2, ncol=17)), # intercepts
#'             cbind(matrix(0, nrow=4, ncol=4), diag(4), -diag(4), matrix(0, nrow=4, ncol=9))) # AR
#'  c <- rep(0, times=6)
#'  Wald_test(fit12, A=A, c=c)
#' @export

Wald_test <- function(stvar, A, c) {
  params <- stvar$params
  stopifnot(is.matrix(A) && ncol(A) == length(params) && nrow(A) <= ncol(A))
  stopifnot(length(c) == nrow(A))
  stopifnot(!is.null(stvar$data))
  if(qr(A)$rank != nrow(A)) stop("The constraint matrix 'A' should have full row rank")

  # Calculate Hessian matrix at the estimate
  Hess <- tryCatch(get_hessian(stvar), error=function(e) {
    message(paste("Failed to calculate Hessian matrix:", e))
    return(NA)})

  # Invert the Hessian matrix
  inv_Hess <- tryCatch(solve(Hess), error=function(e) {
    message(paste("Failed to invert Hessian matrix:", e))
    return(NA)
  })
  if(anyNA(inv_Hess)) stop("Couldn't invert Hessian matrix of the log-likelihood function. Use the function 'get_soc'
                           to make sure the Hessian matrix is positive definite. If it is not pos. def., check the estimates,
                           e.g., with the function 'profile_logliks'.")

  # Calculate the test statistic # t(A%*%params - c)%*%solve(-A%*%inv_Hess%*%t(A))%*%(A%*%params - c)
  test_stat <- as.numeric(crossprod(A%*%params - c, solve(-tcrossprod(A%*%inv_Hess, A), A%*%params - c)))

  # Calculate the p-value
  df <- nrow(A)
  p_value <- pchisq(test_stat, df=df, lower.tail=FALSE)

  # Return
  structure(list(stvar=stvar,
                 A=A,
                 c=c,
                 df=df,
                 test_stat=test_stat,
                 df=df,
                 p_value=p_value,
                 type="Wald test"),
            class="hypotest")
}


#' @title Perform likelihood ratio test for a STVAR model
#'
#' @description \code{LR_test} performs a likelihood ratio test for a STVAR model
#'
#' @param stvar1 an object of class \code{'stvar'} generated by \code{fitSTVAR} or \code{STVAR}, containing
#'   the \strong{freely estimated} model.
#' @param stvar2 an object of class \code{'stvar'} generated by \code{fitSTVAR} or \code{STVAR}, containing
#'   the \strong{constrained} model.
#' @details Performs a likelihood ratio test, testing the null hypothesis that the true parameter value lies
#'   in the constrained parameter space. Under the null, the test statistic is asymptotically
#'   \eqn{\chi^2}-distributed with \eqn{k} degrees of freedom, \eqn{k} being the difference in the dimensions
#'   of the unconstrained and constrained parameter spaces.
#'
#'   \strong{The test is based on the assumption of the standard result of asymptotic normality!}
#'   Also, note that this function does \strong{not} verify that the two models are actually nested.
#' @seealso \code{\link{Wald_test}}, \code{\link{Rao_test}}, \code{\link{fitSTVAR}}, \code{\link{STVAR}},
#'   \code{\link{diagnostic_plot}}, \code{\link{profile_logliks}}, \code{\link{Portmanteau_test}}
#' @inherit Wald_test references return
#' @examples
#' # Logistic Student's t STVAR with p=1, M=2, and the first lag of the second variable
#' # as the switching variable (parameter values were obtained by maximum likelihood estimation;
#'  # fitSTVAR is not used here because the estimation is computationally demanding).
#' params12 <- c(0.62906848, 0.14245295, 2.41245785, 0.66719269, 0.3534745, 0.06041779, -0.34909745,
#'   0.61783824, 0.125769, -0.04094521, -0.99122586, 0.63805416, 0.371575, 0.00314754, 0.03440824,
#'   1.29072533, -0.06067807, 0.18737385, 1.21813844, 5.00884263, 7.70111672)
#' fit12 <- STVAR(data=gdpdef, p=1, M=2, params=params12, weight_function="logistic",
#'   weightfun_pars=c(2, 1), cond_dist="Student")
#' fit12
#'
#' ## Test whether the location parameter equals 1:
#'
#' # Same as the original model but with the location parameter constrained to 1
#' # (parameter values were obtained by maximum likelihood estimation; fitSTVAR
#' # is not used here because the estimation is computationally demanding).
#' params12w <- c(0.6592583, 0.16162866, 1.7811393, 0.38876396, 0.35499367, 0.0576433,
#'  -0.43570508, 0.57337706, 0.16449607, -0.01910167, -0.70747014, 0.75386158, 0.3612087,
#'   0.00241419, 0.03202824, 1.07459924, -0.03432236, 0.14982445, 6.22717097, 8.18575651)
#' fit12w <- STVAR(data=gdpdef, p=1, M=2, params=params12w, weight_function="logistic",
#'  weightfun_pars=c(2, 1), cond_dist="Student",
#'  weight_constraints=list(R=matrix(c(0, 1), nrow=2), r=c(1, 0)))
#'
#' # Test the null hypothesis of the location parameter equal 1:
#' LR_test(fit12, fit12w)
#'
#' ## Test whether the means and AR matrices are identical across the regimes:
#'
#' # Same as the original model but with the mean and AR matrices constrained identical
#' # across the regimes (parameter values were obtained by maximum likelihood estimation;
#' # fitSTVAR is not used here because the estimation is computationally demanding).
#' params12cm <- c(0.76892423, 0.67128089, 0.30824474, 0.03530802, -0.11498402, 0.85942541,
#'  0.39106754, 0.0049437, 0.03897287, 1.44457723, -0.05939876, 0.20885008, 1.23568782,
#'  6.42128475, 7.28733557)
#' fit12cm <- STVAR(data=gdpdef, p=1, M=2, params=params12cm, weight_function="logistic",
#'  weightfun_pars=c(2, 1), parametrization="mean", cond_dist="Student",
#'  mean_constraints=list(1:2), AR_constraints=rbind(diag(4), diag(4)))
#'
#' # Test the null hypothesis of the means and AR matrices being identical across the regimes:
#' LR_test(fit12, fit12cm)
#' @export

LR_test <- function(stvar1, stvar2) {
  stopifnot(class(stvar1) == "stvar")
  stopifnot(class(stvar2) == "stvar")
  stopifnot(length(stvar1$params) > length(stvar2$params))
  stopifnot(stvar1$loglik >= stvar2$loglik)

  test_stat <- as.numeric(2*(stvar1$loglik - stvar2$loglik))
  df <- length(stvar1$params) - length(stvar2$params)
  p_value <- pchisq(test_stat, df=df, lower.tail=FALSE)

  # Return
  structure(list(stvar1=stvar1,
                 stvar2=stvar2,
                 test_stat=test_stat,
                 df=df,
                 p_value=p_value,
                 type="Likelihood ratio test"),
            class="hypotest")
}


#' @title Perform Rao's score test for a STVAR model
#'
#' @description \code{Rao_test} performs Rao's score test for a STVAR model
#'
#' @param stvar an object of class \code{'stvar'} generated by \code{fitSTVAR} or \code{STVAR}, containing
#'   the model specified by the null hypothesis (i.e., \strong{the constrained model}).
#' @details Tests the constraints imposed in the model given in the argument \code{stvar}.
#'   This implementation uses the outer product of gradients approximation in the test statistic.
#'
#'   \strong{The test is based on the assumption of the standard result of asymptotic normality!}
#' @seealso \code{\link{LR_test}}, \code{\link{Wald_test}}, \code{\link{fitSTVAR}}, \code{\link{STVAR}},
#'  \code{\link{diagnostic_plot}}, \code{\link{profile_logliks}}, \code{\link{Portmanteau_test}}
#' @inherit Wald_test references return
#' @examples
#' \donttest{
#' ## These are long running examples that take approximately 10 seconds to run.
#'
#' # Logistic Student's t STVAR with p=1, M=2, and the first lag of the second variable
#' # as the switching variable.
#'
#' ## Test whether the location parameter equal 1:
#'
#' # The model imposing the constraint on the location parameter (parameter values
#' # were obtained by maximum likelihood estimation; fitSTVAR is not used here
#' # because the estimation is computationally demanding):
#' params12w <- c(0.6592583, 0.16162866, 1.7811393, 0.38876396, 0.35499367, 0.0576433,
#'   -0.43570508, 0.57337706, 0.16449607, -0.01910167, -0.70747014, 0.75386158, 0.3612087,
#'   0.00241419, 0.03202824, 1.07459924, -0.03432236, 0.14982445, 6.22717097, 8.18575651)
#' fit12w <- STVAR(data=gdpdef, p=1, M=2, params=params12w, weight_function="logistic",
#'  weightfun_pars=c(2, 1), cond_dist="Student",
#'  weight_constraints=list(R=matrix(c(0, 1), nrow=2), r=c(1, 0)))
#' fit12w
#'
#' # Test the null hypothesis of the location parameter equal 1:
#' Rao_test(fit12w)
#'
#' ## Test whether the means and AR matrices are identical across the regimes:
#'
#' # The model imposing the constraint on the location parameter (parameter values
#' # were obtained by maximum likelihood estimation; fitSTVAR is not used here
#' # because the estimation is computationally demanding):
#' params12cm <- c(0.76892423, 0.67128089, 0.30824474, 0.03530802, -0.11498402, 0.85942541,
#'  0.39106754, 0.0049437, 0.03897287, 1.44457723, -0.05939876, 0.20885008, 1.23568782,
#'  6.42128475, 7.28733557)
#' fit12cm <- STVAR(data=gdpdef, p=1, M=2, params=params12cm, weight_function="logistic",
#'  weightfun_pars=c(2, 1), parametrization="mean", cond_dist="Student",
#'  mean_constraints=list(1:2), AR_constraints=rbind(diag(4), diag(4)))
#'
#' # Test the null hypothesis of the means and AR matrices being identical across the regimes:
#' Rao_test(fit12cm)
#' }
#' @export

Rao_test <- function(stvar) {
  stopifnot(!is.null(stvar$data))

  # Obtain the unconstrained model by expanding the contraints:
  p <- stvar$model$p
  M <- stvar$model$M
  d <- stvar$model$d
  params <- stvar$params
  weight_function <- stvar$model$weight_function
  weightfun_pars <- check_weightfun_pars(data=stvar$data, p=p, M=M, d=d, weight_function=weight_function,
                                         weightfun_pars=stvar$model$weightfun_pars)
  cond_dist <- stvar$model$cond_dist
  parametrization <- stvar$model$parametrization
  identification <- stvar$model$identification
  AR_constraints <- stvar$model$AR_constraints
  mean_constraints <- stvar$model$mean_constraints
  weight_constraints <- stvar$model$weight_constraints
  B_constraints <- stvar$model$B_constraints
  penalized <- stvar$penalized
  penalty_params <- stvar$penalty_params
  allow_unstab <- stvar$allow_unstab
  params <- reform_constrained_pars(p=p, M=M, d=d, params=params,
                                    weight_function=weight_function, weightfun_pars=weightfun_pars,
                                    cond_dist=cond_dist, identification=identification, AR_constraints=AR_constraints,
                                    mean_constraints=mean_constraints, weight_constraints=weight_constraints,
                                    B_constraints=B_constraints)
  new_stvar <- STVAR(data=stvar$data, p=p, M=M, d=d, params=params, parametrization=parametrization,
                     weight_function=weight_function, weightfun_pars=weightfun_pars,
                     cond_dist=cond_dist, identification=identification, AR_constraints=NULL,
                     mean_constraints=NULL, weight_constraints=NULL, B_constraints=NULL,
                     penalized=penalized, penalty_params=penalty_params, allow_unstab=allow_unstab)

  # Calculate the gradient:
  Grad <- tryCatch(get_gradient(new_stvar),
                   error=function(e) {
                     message(paste("Failed to calculate the gradient matrix:", e))
                     return(NA)})
  if(anyNA(Grad)) {
    stop("Couldn't calculate gradient. Check that the estimates are not very close to the boundary of the parameter space.")
  }

  # Gradients of the terms l_t
  foo <- function(x, which_obs=1) {
      # Log-likelihood function as a function of the parameter
      loglikelihood(data=new_stvar$data, p=new_stvar$model$p, M=new_stvar$model$M, params=x,
                    weight_function=new_stvar$model$weight_function,
                    weightfun_pars=new_stvar$model$weightfun_pars,
                    cond_dist=new_stvar$model$cond_dist,
                    parametrization=new_stvar$model$parametrization,
                    identification=new_stvar$model$identification,
                    AR_constraints=new_stvar$model$AR_constraints,
                    mean_constraints=new_stvar$model$mean_constraints,
                    weight_constraints=new_stvar$model$weight_constraints,
                    B_constraints=new_stvar$model$B_constraints,
                    penalized=new_stvar$penalized, penalty_params=new_stvar$penalty_params,
                    allow_unstab=new_stvar$allow_unstab, to_return="terms", minval=NA,)[which_obs]
  }
  T_obs <- nrow(new_stvar$data) - new_stvar$model$p
  outer_prods <- array(dim=c(length(Grad), length(Grad), T_obs))
  # Calculate the outer product of the gradient for each observation, i.e., gradient each term l_t, t=1,...,T:
  for(i1 in 1:T_obs) {
    grad0 <- calc_gradient(x=new_stvar$params, fn=foo, which_obs=i1)
    outer_prods[, , i1] <- tcrossprod(grad0)
  }
  OPG <- apply(outer_prods, MARGIN=1:2, FUN=sum)

  # Invert the OPG matrix
  inv_OPG <- tryCatch(solve(OPG), error=function(e) {
    message(paste("Failed to invert the outer product of gradients matrix:", e))
    return(NA)
  })

 if(anyNA(inv_OPG)) stop("Couldn't invert the outer product of gradients matrix of the log-likelihood function. Make sure the
                          estimates are ok, e.g., with the function 'profile_logliks'.")

  # Calculate the test statistic
  test_stat <- as.numeric(crossprod(Grad, inv_OPG%*%Grad)) # as.numeric(crossprod(Grad, -inv_Hess%*%Grad))

  # Calculate the p-value
  df <- length(new_stvar$params) - length(stvar$params) # The number of constraints
  p_value <- pchisq(test_stat, df=df, lower.tail=FALSE)

  # Return
  structure(list(stvar=stvar,
                 stvar_unconstrained=new_stvar,
                 test_stat=test_stat,
                 df=df,
                 p_value=p_value,
                 type="Rao test"),
            class="hypotest")
}