R/FCVAR_aux.R
In FCVAR: Estimation and Inference for the Fractionally Cointegrated VAR

Documented in FCVARlikeGrid FCVARsim FCVARsimBS FracDiff plot.FCVAR_grid

#' Draw Samples from the FCVAR Model
#'
#' \code{FCVARsim} simulates the FCVAR model as specified by
#' input \code{model} and starting values specified by \code{data}.
#' Errors are drawn from a normal distribution.
#' @param x A \eqn{N x p} matrix of \code{N} starting values for the simulated observations.
#' @param model A list of estimation results, just as if estimated from \code{FCVARest}.
#' The parameters in \code{model} can also be set or adjusted by assigning new values.
#' @param NumPeriods The number of time periods in the simulation.
#' @return A \code{NumPeriods} by \eqn{p} matrix \code{xBS} of simulated observations.
#' @examples
#' \donttest{
#' opt <- FCVARoptions()
#' opt$gridSearch   <- 0 # Disable grid search in optimization.
#' opt$dbMin        <- c(0.01, 0.01) # Set lower bound for d,b.
#' opt$dbMax        <- c(2.00, 2.00) # Set upper bound for d,b.
#' opt$constrained  <- 0 # Impose restriction dbMax >= d >= b >= dbMin ? 1 <- yes, 0 <- no.
#' x <- votingJNP2014[, c("lib", "ir_can", "un_can")]
#' results <- FCVARestn(x, k = 2, r = 1, opt)
#' x_sim <- FCVARsim(x[1:10, ], results, NumPeriods = 100)
#' }
#' @family FCVAR auxiliary functions
#' @seealso \code{FCVARoptions} to set default estimation options.
#' \code{FCVARestn} for the specification of the \code{model}.
#' Use \code{FCVARsim} to draw a sample from the FCVAR model.
#' For simulations intended for bootstrapping statistics, use \code{FCVARsimBS}.
#' @export
#'
FCVARsim <- function(x, model, NumPeriods) {


  #--------------------------------------------------------------------------------
  # Preliminary definitions
  #--------------------------------------------------------------------------------

  p <- ncol(x)
  opt <- model$options
  cf  <- model$coeffs
  d <- cf$db[1]
  b <- cf$db[2]

  # Generate disturbance term
  err <- matrix(stats::rnorm(NumPeriods*p), nrow = NumPeriods, ncol = p)

  #--------------------------------------------------------------------------------
  # Recursively generate simulated data values
  #--------------------------------------------------------------------------------

  # Initialize simulation with starting values.
  x_sim <- x
  for (i in 1:NumPeriods) {


    # Append x_sim with zeros to simplify calculations.
    x_sim <- rbind(x_sim, rep(0, p))
    T <- nrow(x_sim)

    # Adjust by level parameter if present.
    if(opt$levelParam) {
      y <- x_sim - matrix(1, nrow = T, ncol = 1) %*% cf$muHat
    } else {
      y <- x_sim
    }


    # Main term, take fractional lag.
    z <- Lbk(y,d,1)

    # Error correction term.
    if (!is.null(cf$alphaHat)) {

      z <- z + FracDiff( Lbk(y, b, 1), d - b ) %*% t(cf$PiHat)

      if (opt$rConstant) {
        z <- z + FracDiff( Lbk(matrix(1, nrow = T, ncol = 1), b, 1), d - b ) %*%
          cf$rhoHat %*% t(cf$alphaHat)
      }

    }


    # Add unrestricted constant if present.
    if(opt$unrConstant) {
      z <- z + matrix(1, nrow = T, ncol = 1) %*% t(cf$xiHat)
    }


    # Add lags if present.
    if (!is.null(cf$GammaHat)) {
      k <- ncol(cf$GammaHat)/p
      z <- z +  FracDiff(  Lbk( y , b, k)  , d) %*% t(cf$GammaHat)
    }


    # Adjust by level parameter if present.
    if (opt$levelParam) {
      z <- z + matrix(1, nrow = T, ncol = 1) %*% cf$muHat
    }


    # Add disturbance term
    z[T,] <- z[T,] + err[i, ]

    # Append to x_sim matrix.
    x_sim <- rbind(x_sim[1:(T-1), ], z[T,])

  }

  #--------------------------------------------------------------------------------
  # Return simulated data values, including (EXcluding?) initial values.
  #--------------------------------------------------------------------------------

  x_sim <- x_sim[(nrow(x)+1):nrow(x_sim), ]

  return(x_sim)
}


#' Draw Bootstrap Samples from the FCVAR Model
#'
#' \code{FCVARsimBS} simulates the FCVAR model as specified by
#' input \code{model} and starting values specified by \code{data}.
#' It creates a wild bootstrap sample by augmenting each iteration
#' with a bootstrap error. The errors are sampled from the
#' residuals specified under the \code{model} input and have a
#' positive or negative sign with equal probability (the Rademacher distribution).
#' @param data A \eqn{T x p} matrix of starting values for the simulated realizations.
#' @param model A list of estimation results, just as if estimated from \code{FCVARest}.
#' The parameters in \code{model} can also be set or adjusted by assigning new values.
#' @param NumPeriods The number of time periods in the simulation.
#' @return A \code{NumPeriods} by \eqn{p} matrix \code{xBS} of simulated bootstrap values.
#' @examples
#' \donttest{
#' opt <- FCVARoptions()
#' opt$gridSearch   <- 0 # Disable grid search in optimization.
#' opt$dbMin        <- c(0.01, 0.01) # Set lower bound for d,b.
#' opt$dbMax        <- c(2.00, 2.00) # Set upper bound for d,b.
#' opt$constrained  <- 0 # Impose restriction dbMax >= d >= b >= dbMin ? 1 <- yes, 0 <- no.
#' x <- votingJNP2014[, c("lib", "ir_can", "un_can")]
#' results <- FCVARestn(x, k = 2, r = 1, opt)
#' xBS <- FCVARsimBS(x[1:10, ], results, NumPeriods = 100)
#' }
#' @family FCVAR auxiliary functions
#' @seealso \code{FCVARoptions} to set default estimation options.
#' \code{FCVARestn} for the specification of the \code{model}.
#' Use \code{FCVARsim} to draw a sample from the FCVAR model.
#' For simulations intended for bootstrapping statistics, use \code{FCVARsimBS}.
#' @export
#'
FCVARsimBS <- function(data, model, NumPeriods) {


  #--------------------------------------------------------------------------------
  # Preliminary definitions
  #--------------------------------------------------------------------------------

  x <- data
  p <- ncol(data)
  opt <- model$options
  cf  <- model$coeffs
  d <- cf$db[1]
  b <- cf$db[2]
  T <- nrow(model$Residuals)

  # Generate disturbance term for use in the bootstrap.

  # Center residuals
  res <- model$Residuals -
    matrix(1, nrow = nrow(model$Residuals), ncol = 1) %*% colMeans(model$Residuals)

  # Generate draws from Rademacher distribution for Wild bootstrap
  eRD <- - matrix(1, nrow = T, ncol = p) +
    2*( matrix(stats::rnorm(T), nrow = T, ncol = 1) > 0) %*% matrix(1, nrow = 1, ncol = p)

  # Generate error term
  err <- res * eRD



  #--------------------------------------------------------------------------------
  # Recursively generate bootstrap sample
  #--------------------------------------------------------------------------------


  for (i in 1:NumPeriods) {


    # Append x with zeros to simplify calculations
    x <- rbind(x, rep(0, p))
    T <- nrow(x)

    # Adjust by level parameter if present
    if(opt$levelParam) {
      y <- x - matrix(1, nrow = T, ncol = 1) %*% cf$muHat
    } else {
      y <- x
    }


    # Main term, take fractional lag
    z <- Lbk(y, d, 1)


    # Error correction term
    if( !is.null(cf$alphaHat) && !is.na(cf$alphaHat[1])) {
      z <- z + FracDiff( Lbk(y, b, 1), d - b ) %*% t(cf$PiHat)
      if(opt$rConstant) {
        z <- z + FracDiff( Lbk(matrix(1, nrow = T, ncol = 1), b, 1), d - b ) %*%
          cf$rhoHat %*% t(cf$alphaHat)
      }

    }


    # Add unrestricted constant if present
    if(opt$unrConstant) {
      z <- z + matrix(1, nrow = T, ncol = 1) %*% t(cf$xiHat)
    }


    # Add lags if present
    if(!is.null(cf$GammaHat)) {
      k <- ncol(cf$GammaHat)/p
      z <- z +  FracDiff(  Lbk( y , b, k)  , d) %*% t(cf$GammaHat)
    }


    # Adjust by level parameter if present
    if(opt$levelParam) {
      z <- z + matrix(1, nrow = T,ncol = 1) %*% cf$muHat
    }


    # Add disturbance term
    z[T, ] <- z[T, ] + err[i, ]

    # Append generated observation to x matrix
    x <- rbind(x[1:(T-1), ], z[T, ])

  }


  #--------------------------------------------------------------------------------
  # Return bootstrap sample
  #--------------------------------------------------------------------------------

  xBS <- x[(nrow(data)+1):nrow(x), ]

  return(xBS)
}


# Calculate Parameters of the Cointegrated VAR Model
#
# \code{GetParams} calculates the parameters of the cointegrated VAR model
#  from a multivariate sample.
#  This function uses FWL and reduced rank regression to obtain
# the estimates of \code{Alpha}, \code{Beta}, \code{Rho}, \code{Pi},
# \code{Gamma}, and \code{Omega}, using the notation in Johansen and Nielsen (2012).
# Note that Roxygen comments are excluded to keep this function internal.
# However, the contents of Roxygen comments are shown below for those who read the scripts.
#
# @param x A matrix of variables to be included in the system.
# @param k The number of lags in the system.
# @param r The cointegrating rank.
# @param db The orders of fractional integration.
# @param opt An S3 object of class \code{FCVAR_opt} that stores the chosen estimation options,
# generated from \code{FCVARoptions()}.
# @return A list object \code{estimates} containing the estimates,
# including the following parameters:
# \describe{
#   \item{\code{db}}{The orders of fractional integration (taken directly from the input).}
#   \item{\code{alphaHat}}{A \eqn{p x r} matrix of adjustment parameters.}
#   \item{\code{betaHat}}{A \eqn{p x r} matrix of cointegrating vectors.
#       The \eqn{r x 1} vector \eqn{\beta x_t} is the stationary cointegration relations.}
#   \item{\code{rhoHat}}{A \eqn{p x 1} vector of restricted constatnts.}
#   \item{\code{piHat}}{A \eqn{p x p} matrix \eqn{\Pi = \alpha \beta'} of long-run levels.}
#   \item{\code{OmegaHat}}{A \eqn{p x p} covariance matrix of the error terms.}
#   \item{\code{GammaHat}}{A ( \eqn{p x kp} matrix \code{cbind(GammaHat1,...,GammaHatk)})
#   of autoregressive coefficients. }
# }
# @examples
# opt <- FCVARoptions()
# x <- votingJNP2014[, c("lib", "ir_can", "un_can")]
# estimates <- GetParams(x, k = 2, r = 1, db = c(1, 1), opt)
# @family FCVAR auxiliary functions
# @seealso \code{FCVARoptions} to set default estimation options.
# \code{FCVARestn} calls \code{GetParams} to estimate the FCVAR model.
# @references Johansen, S. and M. \enc{Ø}{O}. Nielsen (2012).
# "Likelihood inference for a fractionally cointegrated
# vector autoregressive model," Econometrica 80, 2667-2732.
# @export
#
GetParams <- function(x, k, r, db, opt) {


  Z_array <- TransformData(x, k, db, opt)
  Z0 <- Z_array$Z0
  Z1 <- Z_array$Z1
  Z2 <- Z_array$Z2
  Z3 <- Z_array$Z3

  T <- nrow(Z0)
  p <- ncol(Z0)
  p1 <- ncol(Z1)


  #--------------------------------------------------------------------------------
  # Concentrate out the unrestricted constant if present
  #--------------------------------------------------------------------------------

  if(opt$unrConstant) {
    # Note Z3*inv(Z3*Z3)*Z3 is just a matrix of ones / T
    Z0hat <- Z0 - Z3 %*% ( solve(t(Z3) %*% Z3) %*% t(Z3) %*% Z0 )
    Z1hat <- Z1 - Z3 %*% ( solve(t(Z3) %*% Z3) %*% t(Z3) %*% Z1 )
    if(k > 0) {
      Z2hat <- Z2 - Z3 %*% ( solve(t(Z3) %*% Z3) %*% t(Z3) %*% Z2 )
    } else {
      Z2hat <- Z2
    }

  } else {
    Z0hat <- Z0
    Z1hat <- Z1
    Z2hat <- Z2
  }

  #--------------------------------------------------------------------------------
  # FWL Regressions
  #--------------------------------------------------------------------------------

  if ( (k == 0) || is.infinite(kappa(t(Z2hat) %*% Z2hat)) ) {

    # No lags, so there are no effects of Z2.
    R0 <- Z0hat
    R1 <- Z1hat
    # Includes the case when Z2hat is near zero and
    # t(Z2hat) %*% Z2hat is ill-conditioned.

  } else {

    # Lags included: Obtain the residuals from FWL regressions.
    # Unless, of course, Z2hat is near zero and
    # t(Z2hat) %*% Z2hat is ill-conditioned.
    R0 <- Z0hat - Z2hat %*% ( solve(t(Z2hat) %*% Z2hat) %*% t(Z2hat) %*% Z0hat )
    R1 <- Z1hat - Z2hat %*% ( solve(t(Z2hat) %*% Z2hat) %*% t(Z2hat) %*% Z1hat )
  }


  # Calculate Sij matrices for reduced rank regression.
  S00 <- t(R0) %*% R0/T
  S01 <- t(R0) %*% R1/T
  S10 <- t(R1) %*% R0/T
  S11 <- t(R1) %*% R1/T

  # Calculate reduced rank estimate of Pi.
  if (r == 0) {

    betaHat <- NULL
    betaStar <- NA
    alphaHat <- NA
    PiHat <- NULL
    rhoHat <- NULL
    OmegaHat <- S00

  } else if (( r > 0 ) & ( r < p )) {


    # Calculate solution from eigenvectors
    eig_out <- eigen( solve(S11) %*% S10 %*% solve(S00) %*% S01 )
    D <- eig_out$values # Only the vector, not the diagonal matrix.
    V1 <- eig_out$vectors

    V2 <- t( cbind(t(V1), D)[order(D), ] )

    V <- V2[1:p1, , drop = FALSE]
    betaStar <- V[ 1:p1, seq(p1, p1-r+1, by = -1), drop = FALSE]



    # If either alpha or beta is restricted, then the likelihood is
    #   maximized subject to the constraints. This section of the code
    #   replaces the call to the switching algorithm in the previous
    #   version.
    if (!is.null(opt$R_Alpha) | !is.null(opt$R_Beta) ) {


      switched_mats <- GetRestrictedParams(betaStar, S00, S01, S11, T, p, opt)
      betaStar <- switched_mats$betaStar
      alphaHat <- switched_mats$alphaHat
      OmegaHat <- switched_mats$OmegaHat

      betaHat <- betaStar[1:p, 1:r, drop = FALSE]
      PiHat <- alphaHat %*% t(betaHat)



    } else {

      # Otherwise, alpha and beta are unrestricted, but unidentified.
      alphaHat <- S01 %*% betaStar %*% solve(t(betaStar) %*% S11 %*% betaStar)
      OmegaHat <- S00 - alphaHat %*% t(betaStar) %*% S11 %*% betaStar %*% t(alphaHat)
      betaHat <- matrix(betaStar[1:p, 1:r], nrow = p, ncol = r)
      PiHat <- alphaHat %*% t(betaHat)
      # Transform betaHat and alphaHat to identify beta.
      #   The G matrix is used to identify beta by filling the first
      #   rxr block of betaHat with an identity matrix.
      G <- solve(betaHat[1:r,1:r])
      betaHat <- betaHat %*% G
      betaStar <- betaStar %*% G
      # alphaHat is post multiplied by G^-1 so that Pi= a(G^-1)Gb <- ab
      alphaHat <- alphaHat %*% t(solve(G))

    }



    # Extract coefficient vector for restricted constant model if required.
    if (opt$rConstant) {
      rhoHat <- betaStar[p1, ]
    } else {
      rhoHat <- NULL
    }


  } else {

    # (r = p) and do no need reduced rank regression
    V <- S01 %*% solve(S11)
    betaHat <- t(V[, 1:p, drop = FALSE])
    # For restriced constant, rho <- last column of V.
    alphaHat <- diag(p)
    PiHat <- betaHat
    OmegaHat <- S00 - S01 %*% solve(S11) %*% S10

    # Extract coefficient vector for restricted constant model if required.
    if (opt$rConstant) {
      rhoHat <- t(V[ , p1])
    } else {
      rhoHat <- NULL
    }

    # Check conformability:
    betaStar <- rbind(betaHat, rhoHat)

  }

  # Calculate PiStar independently of how betaStar was estimated.
  PiStar <- alphaHat %*% t(betaStar)

  # Perform OLS regressions to calculate unrestricted constant and
  #   GammaHat matrices if necessary.
  xiHat <- NULL
  if (k == 0) {
    GammaHat <- NULL
  } else {

    if(r > 0) {
      GammaHat <- t( solve(t(Z2hat) %*% Z2hat) %*% t(Z2hat) %*% (Z0hat - Z1hat %*% t(PiStar)) )
    } else {
      GammaHat <- t( solve(t(Z2hat) %*% Z2hat) %*% t(Z2hat) %*% Z0hat )
    }

  }



  if(opt$unrConstant==1) {
    if(k>0) {
      if(r>0) {
        xiHat <- t( solve(t(Z3) %*% Z3) %*% t(Z3) %*% (Z0 - Z1 %*% t(PiStar) - Z2 %*% t(GammaHat)) )
      } else {
        xiHat <- t( solve(t(Z3) %*% Z3) %*% t(Z3) %*% (Z0 - Z2 %*% t(GammaHat)) )
      }

    } else {
      if(r>0) {
        xiHat <- t( solve(t(Z3) %*% Z3) %*% t(Z3) %*% (Z0 - Z1 %*% t(PiStar)) )
      } else {
        xiHat <- t( solve(t(Z3) %*% Z3) %*% t(Z3) %*% Z0 )
      }

    }

  }




  #--------------------------------------------------------------------------------
  # Return the results in a list
  #--------------------------------------------------------------------------------

  estimates <- list(

    db = db,
    alphaHat = alphaHat,
    betaHat = betaHat,
    rhoHat = rhoHat,
    PiHat = PiHat,
    OmegaHat = OmegaHat,
    GammaHat = GammaHat,
    xiHat = xiHat

  )


  return(estimates)

}


#' Grid Search to Maximize Likelihood Function
#'
#' \code{FCVARlikeGrid} performs a grid-search optimization
#' by calculating the likelihood function
#' on a grid of candidate parameter values.
#' This function evaluates the likelihood over a grid of values
#' 	for \code{c(d,b)} (or \code{phi}).
#' 	It can be used when parameter estimates are sensitive to
#' 	starting values to give an approximation of the global maximum that can
#' 	then be used as the starting value in the numerical optimization in
#' 	\code{FCVARestn}.
#' 	\code{plot.FCVAR_grid} plots the likelihood function from \code{FCVARlikeGrid}.
#'
#' @param x A matrix of variables to be included in the system.
#' @param k The number of lags in the system.
#' @param r The cointegrating rank.
#' @param opt An S3 object of class \code{FCVAR_opt} that stores the chosen estimation options,
#' generated from \code{FCVARoptions()}.
#' @return An S3 object of type \code{FCVAR_grid} containing the optimization results,
#' including the following parameters:
#' \describe{
#'   \item{\code{params}}{A vector \code{params} of \code{d} and \code{b}
#'   (and \code{mu} if level parameter is selected)
#'   corresponding to a maximum over the grid of \code{c(d,b)} or \code{phi}.}
#'   \item{\code{dbHatStar}}{A vector of \code{d} and \code{b}
#'   corresponding to a maximum over the grid of \code{c(d,b)} or \code{phi}.}
#'   \item{\code{muHatStar}}{A vector of the optimal \code{mu} if level parameter is selected. }
#'   \item{\code{Grid2d}}{An indicator for whether or not the optimization
#'   is conducted over a 2-dimensional parameter space,
#'   i.e. if there is no equality restriction on \code{d} and \code{b}.}
#'   \item{\code{dGrid}}{A vector of the grid points in the parameter \code{d},
#'     after any transformations for restrictions, if any.}
#'   \item{\code{bGrid}}{A vector of the grid points in the parameter \code{b},
#'     after any transformations for restrictions, if any.}
#'   \item{\code{dGrid_orig}}{A vector of the grid points in the parameter \code{d},
#'     in units of the fractional integration parameter.}
#'   \item{\code{bGrid_orig}}{A vector of the grid points in the parameter \code{b},
#'     in units of the fractional integration parameter.}
#'   \item{\code{like}}{The maximum value of the likelihood function over the chosen grid.}
#'   \item{\code{k}}{The number of lags in the system.}
#'   \item{\code{r}}{The cointegrating rank.}
#'   \item{\code{opt}}{An S3 object of class \code{FCVAR_opt} that stores the chosen estimation options,
#'     generated from \code{FCVARoptions()}.}
#' }
#'
#' @examples
#' # Restrict equality of fractional parameters.
#' \donttest{
#' opt <- FCVARoptions()
#' opt$dbStep1D     <- 0.2 # Coarser grid for plotting example.
#' opt$dbMin        <- c(0.01, 0.01) # Set lower bound for d,b.
#' opt$dbMax        <- c(2.00, 2.00) # Set upper bound for d,b.
#' opt$constrained  <- 0 # impose restriction dbMax >= d >= b >= dbMin ? 1 <- yes, 0 <- no.
#' opt$restrictDB   <- 1 # impose restriction d=b ? 1 <- yes, 0 <- no.
#' opt$progress     <- 2 # Show progress report on each value of b.
#' x <- votingJNP2014[, c("lib", "ir_can", "un_can")]
#' likeGrid_params <- FCVARlikeGrid(x, k = 2, r = 1, opt)
#' plot(likeGrid_params)
#' }
#'
#' # Linear restriction on fractional parameters.
#' \donttest{
#' opt <- FCVARoptions()
#' opt$dbStep1D     <- 0.2 # Coarser grid for plotting example.
#' opt$dbMin        <- c(0.01, 0.01) # Set lower bound for d,b.
#' opt$dbMax        <- c(2.00, 2.00) # Set upper bound for d,b.
#' opt$constrained  <- 0 # impose restriction dbMax >= d >= b >= dbMin ? 1 <- yes, 0 <- no.
#' opt$restrictDB   <- 0 # impose restriction d=b ? 1 <- yes, 0 <- no.
#' # Impose linear restriction on d and b:
#' opt$R_psi        <- matrix(c(2, -1), nrow = 1, ncol = 2)
#' opt$r_psi        <- 0.5
#' opt$progress     <- 2 # Show progress report on each value of b.
#' x <- votingJNP2014[, c("lib", "ir_can", "un_can")]
#' likeGrid_params <- FCVARlikeGrid(x, k = 2, r = 1, opt)
#' plot(likeGrid_params)
#' }
#'
#' # Constrained 2-dimensional optimization.
#' # Impose restriction dbMax >= d >= b >= dbMin.
#' \donttest{
#' opt <- FCVARoptions()
#' opt$dbStep1D     <- 0.2 # Coarser grid for plotting example.
#' opt$dbStep2D     <- 0.2 # Coarser grid for plotting example.
#' opt$dbMin        <- c(0.01, 0.01) # Set lower bound for d,b.
#' opt$dbMax        <- c(2.00, 2.00) # Set upper bound for d,b.
#' opt$constrained  <- 1 # impose restriction dbMax >= d >= b >= dbMin ? 1 <- yes, 0 <- no.
#' opt$restrictDB   <- 0 # impose restriction d=b ? 1 <- yes, 0 <- no.
#' opt$progress     <- 2 # Show progress report on each value of b.
#' x <- votingJNP2014[, c("lib", "ir_can", "un_can")]
#' likeGrid_params <- FCVARlikeGrid(x, k = 2, r = 1, opt)
#' }
#'
#' # Unconstrained 2-dimensional optimization.
#' \donttest{
#' opt <- FCVARoptions()
#' opt$dbStep1D     <- 0.1 # Coarser grid for plotting example.
#' opt$dbStep2D     <- 0.2 # Coarser grid for plotting example.
#' opt$dbMin        <- c(0.01, 0.01) # Set lower bound for d,b.
#' opt$dbMax        <- c(2.00, 2.00) # Set upper bound for d,b.
#' opt$constrained  <- 0 # impose restriction dbMax >= d >= b >= dbMin ? 1 <- yes, 0 <- no.
#' opt$restrictDB   <- 0 # impose restriction d=b ? 1 <- yes, 0 <- no.
#' opt$progress     <- 2 # Show progress report on each value of b.
#' x <- votingJNP2014[, c("lib", "ir_can", "un_can")]
#' likeGrid_params <- FCVARlikeGrid(x, k = 2, r = 1, opt)
#' }
#' @family FCVAR auxiliary functions
#' @seealso \code{FCVARoptions} to set default estimation options.
#' \code{plot.FCVAR_grid} plots the likelihood function from \code{FCVARlikeGrid}.
#' @note If \code{opt$LocalMax == 0}, \code{FCVARlikeGrid} returns the parameter values
#'       corresponding to the global maximum of the likelihood on the grid.
#'       If \code{opt$LocalMax == 1}, \code{FCVARlikeGrid} returns the parameter values for the
#'       local maximum corresponding to the highest value of \code{b}. This
#'       alleviates the identification problem mentioned in Johansen and
#'       Nielsen (2010, section 2.3).
#' @references Johansen, S. and M. \enc{Ø}{O}. Nielsen (2010).
#' "Likelihood inference for a nonstationary fractional
#' autoregressive model," Journal of Econometrics 158, 51-66.
#' @export
#'
FCVARlikeGrid <- function(x, k, r, opt) {

  p <- ncol(x)

  # Update options to account for all settings and restrictions.
  # And to output warning or error messages if opt is incorrectly specified.
  opt <- FCVARoptionUpdates(opt, p, r)


  #--------------------------------------------------------------------------------
  # INITIALIZATION
  #--------------------------------------------------------------------------------

  # Check if performing a 2-dimensional or 1-dimensional grid search. The
  #   step size will be smaller if searching in only one dimension.
  if(is.null(opt$R_psi)) {
    Grid2d <- 1
    dbStep <- opt$dbStep2D # This variable was added to options.
  } else {
    Grid2d <- 0
    dbStep <- opt$dbStep1D
  }



  # If equality restrictions are imposed, need to construct a grid over
  #   phi and obtain db <- H*phi + h.
  if(!is.null(opt$R_psi)) {
    # This set of variables is defined for easy translation between
    #   phi (unrestricted parameter) and d,b (restricted parameters).
    R <- opt$R_psi
    s <- opt$r_psi
    H <- pracma::null(R)
    h <- t(R) %*% solve(R %*% t(R)) %*% s
  }


  # GetBounds returns upper and lower bounds in terms of phi when
  # restrictions are imposed and in terms of d and b when they are
  # unrestricted. It also adjusts for limits on d and b imposed by user.
  UB_LB_bounds <- GetBounds(opt)
  dbMax <- UB_LB_bounds$UB
  dbMin <- UB_LB_bounds$LB



  # Set up the grid.
  if(Grid2d) {

    # Search over d as well as b.
    # Set up grid for d parameter
    dGrid <-  seq(dbMin[1], dbMax[1], by = dbStep)
    # Number of grid points
    nD    <- length(dGrid)
    # Set up grid for b parameter
    bGrid <-  seq(dbMin[2], dbMax[2], by = dbStep)
    # Number of grid points
    nB    <- length(bGrid)
    if(opt$constrained) {
      # Since it is possible to have different grid lengths for d and
      # b when the parameters have different max/min values, the
      # following block calculates the total number of iterations
      # by counting the number of instances in which the grid values
      # in b are less than or equal to those in d. This should be
      # moved to the beginning of the loop for faster computation
      # later.
      bdGrid <- expand.grid(dGrid = dGrid, bGrid = bGrid)
      totIters <- sum(bdGrid[, 'bGrid'] <= bdGrid[, 'dGrid'])
    } else {
      # Unconstrained so search through entire grid for d.
      dStart   <- 1
      totIters <- nB*nD
    }


  } else {
    # Only searching over one parameter.
    bGrid <-  seq(dbMin[1], dbMax[1], by = dbStep)
    dGrid <-  NA
    nB    <- length(bGrid)
    nD    <- 1
    dStart   <- 1
    totIters <- nB
  }

  # Initialize progress reports.
  if(opt$progress != 0) {
    if(opt$progress == 1) {
      progbar <- utils::txtProgressBar(min = 0, max = totIters, style = 3)

    }

  }


  # Create a matrix of NA's to store likelihoods, we use NA's here
  #   because NA entries are not plotted and do not affect the finding
  #   of the maximum.
  like  <- matrix(NA, nrow = nB, ncol = nD)

  # Initialize storage bin and starting values for optimization involving
  #   level parameter.
  if(opt$levelParam) {
    mu  <- array(0, dim = c(p, nB, nD))
    StartVal <- x[1, ,drop = FALSE]
  }


  #--------------------------------------------------------------------------------
  # CALCULATE LIKELIHOOD OVER THE GRID
  #--------------------------------------------------------------------------------



  iterCount <- 0

  for (iB in 1:nB) {

    b <- bGrid[iB]

    # Print a progress report of the preferred style.
    if (opt$progress == 1) {
      utils::setTxtProgressBar(progbar, value = iterCount)
    } else if (opt$progress == 2) {
      message(sprintf('Now estimating for iteration %d of %d: b = %f.',
                      iterCount, totIters, b))
    }

    # If d>=b (constrained) then search only over values of d that are
    #   greater than or equal to b. Also, perform this operation only if
    #   searching over both d and b.
    if(opt$constrained & Grid2d) {
      # Start at the index that corresponds to the first value in the
      # grid for d that is >= b
      dStart <- which(dGrid >= b)[1]
    }


    for (iD in dStart:nD) {

      iterCount <- iterCount + 1

      # db is defined in terms of d,b for use by FCVARlikeMU
      # if level parameters are present and for displaying in
      # the output. phi is used by FCVARlike, which can
      # handle both phi or d,b and makes appropriate
      # adjustments inside the function.
      if(is.null(opt$R_psi)) {
        d <- dGrid[iD]
        db <- cbind(d, b)
        phi <- db
      } else {
        phi <- bGrid[iB]
        db <- H*phi + h
      }


      if(opt$levelParam) {
        # Optimize over level parameter for given (d,b).
        # [ muHat, maxLike, ! ] ...
        # <- fminunc(@( params ) -FCVARlikeMu(x, db, params, k, r, opt), ...
        #            StartVal, opt$UncFminOptions )

        min_out <- stats::optim(StartVal,
                         function(params) {-FCVARlikeMu(params, x, db, k, r, opt)})

        muHat <- min_out$par
        maxLike <- min_out$value

        # Store the results.
        like[iB,iD] <- -maxLike
        mu[, iB,iD] <- muHat

      } else {
        # Only called if no level parameters. phi contains
        # either one or two parameters, depending on
        # whether or not restrictions have been imposed.
        like[iB,iD] <- FCVARlike(phi, x, k, r, opt)
      }


    }


  }


  #--------------------------------------------------------------------------------
  # FIND THE MAX OVER THE GRID
  #--------------------------------------------------------------------------------

  if(opt$LocalMax) {

    # Local max.
    if(Grid2d) {

      # With constrained model, with d >= b, need to pad like with low values
      # where b > d, where the like would otherwise be NA,
      # since they weren't calculated.
      if(opt$constrained) {
        like_padded <- like
        min_like_out_bds = min(like, na.rm = TRUE) - 1000
        for (iB in 2:nB) {
          # END at the index that corresponds to the LAST value in the
          # grid for d that is < b.
          b <- bGrid[iB]
          d_lt_b_grid <- which(dGrid < b)
          dEnd <- d_lt_b_grid[length(d_lt_b_grid)]
          for (iD in 1:dEnd) {
            like_padded[iB,iD] <- min_like_out_bds - bGrid[iB]*100 + dGrid[iD]*100
          }
        }
      } else {
        like_padded <- like
      }


      # Implement local max:
      loc_max_out <- find_local_max(as.array(like_padded))
      indexB <- loc_max_out$row
      indexD <- loc_max_out$col

      # Store value of local max.
      local_max <- list(b = bGrid[indexB],
                      d = dGrid[indexD],
                      like = loc_max_out$value)

      # Record temporary value of maximal likelihood value.
      # It may be replaced to avoid the identification problem.
      max_like <- max(like, na.rm = TRUE)

    } else {

      # Implement local max:
      # Need to pad 1-D array with boundary columns.
      loc_max_out <- find_local_max(as.array(cbind(like - 1,
                                                   like,
                                                   like - 1)))
      indexB <- loc_max_out$row
      indexD <- 1

      # Store value of local max.
      local_max <- list(b = bGrid[indexB],
                        d = rep(NA, length(bGrid[indexB])),
                        like = loc_max_out$value)

      # Record temporary value of maximal likelihood value.
      # It may be replaced to avoid the identification problem.
      max_like <- max(like, na.rm = TRUE)
    }


    # If there is no local max, return global max.
    if(is.null(indexB) | is.null(indexD)) {
      warning('Failure to find local maximum. Returning global maximum instead.',
              'Check whether global maximum is on a boundary and, if so, ',
              'consider expanding the constraints on d and b. ')

      max_like <- max(like, na.rm = TRUE)
      indexB_and_D <- which(like == max_like, arr.ind = TRUE)

      indexB <- indexB_and_D[, 1, drop = FALSE]
      indexD <- indexB_and_D[, 2, drop = FALSE]
    }


  } else {
    # Global max.
    max_like <- max(like, na.rm = TRUE)
    indexB_and_D <- which(like == max(like, na.rm = TRUE), arr.ind = TRUE)

    indexB <- indexB_and_D[, 1, drop = FALSE]
    indexD <- indexB_and_D[, 2, drop = FALSE]
  }


  # Choose value among maxima, if not unique.
  if(length(indexD) > 1 | length(indexB) > 1) {
    # If maximum is not unique, take the index pair corresponding to
    # the highest value of b.
    if(Grid2d) {
      # Sort in ascending order according to indexB.
      # In a 2-D grid, the ordering is unambiguous:
      # sorting in increasing indexB is also increasing in b.
      indBindx <- order(indexB)
      indexB <- indexB[indBindx]
      indexD <- indexD[indBindx]

      # With 2-D grid, there is no transformation,
      # so b is monotinically related to indexB.
      sign_H2 <- 1

    } else {
      # Sort in ascending order by b.
      # indexB <- indexB[order(indexB)]

      # Not so fast!
      # In a 1-D grid, the ordering is also unambiguous
      # WHEN THERE ARE NO RESTRICTIONS, aside from d = b:
      # sorting in increasing indexB is also increasing in b.
      # However, when a linear restriction is placed on d and b,
      # the grid is on the 1-D parameter phi,
      # and the relationship with b is either positive or negative,
      # depending on the second element of H in H*[d, b]' + h.
      if(!is.null(opt$R_psi)) {
        sign_H2 <- sign(H[2])
      } else {
        sign_H2 <- 1
      }
      # Then, sort in the appropriate order.
      indexB <- indexB[order(sign_H2*indexB)]

    }

    # Take last value, which is the local maximum with highest value of b.
    indexB <- indexB[length(indexB)]
    indexD <- indexD[length(indexD)] # Works even in 1-D grid, since both ordered by b.

    # Choosing this local maximum avoids the identification problem.
    # max_like <- local_max$like[which.max(local_max$b)]
    # Again, the order depends on the sign of H[2], if it is a restricted model.
    max_like <- local_max$like[which.max(sign_H2*local_max$b)]

    # cat(sprintf('\nWarning, grid search did not find a unique maximum of the log-likelihood function.\n')  )
    message('Grid search did not find a unique maximum of the log-likelihood function.',
            '\nInspect the plot of the log-likelihood function to verify your choice of the LocalMax option.')
  }

  # Translate to d,b.
  if(!is.null(opt$R_psi)) {
    # Translate the parameter estimates.
    dbHatStar <- t(H %*% bGrid[indexB] + h)

    # Translate the grid points to units of the original parameters.
    # This is useful for plotting the likelihood function.
    bGrid_orig <- 0*bGrid
    dGrid_orig <- 0*dGrid
    for (i in 1:length(bGrid)) {
      dbHatStar_i <- t(H %*% bGrid[i] + h)
      dGrid_orig[i] <- dbHatStar_i[1]
      bGrid_orig[i] <- dbHatStar_i[2]
    }

    # Translate local maxima.
    for (i in 1:length(local_max$b)) {
      # Recall that, in this case, bGrid is actually the grid for phi,
      # the 1-D search parameter in the restricted definition of [d, b].
      local_max_db <- matrix(c(local_max$b[i]), nrow = length(local_max$b[i]), ncol = 1)

      dbHatStar_i <- t(H %*% local_max_db + h)
      # These were switched:
      # local_max$b[i] <- dbHatStar_i[1]
      # local_max$d[i] <- dbHatStar_i[2]
      # d is the first parameter in [d, b]; b is the second.
      local_max$d[i] <- dbHatStar_i[1]
      local_max$b[i] <- dbHatStar_i[2]
    }

  } else {
    # No translation necessary.
    dbHatStar <- cbind(dGrid[indexD], bGrid[indexB])
    bGrid_orig <- bGrid
    dGrid_orig <- dGrid
  }



  # Print final progress report.
  if (opt$progress != 0) {
    if ( iterCount == totIters) {
      if (opt$progress == 1) {
        utils::setTxtProgressBar(progbar, value = iterCount)
        # message(sprintf('\nProgress : %5.1f%%, b = %4.2f, d = %4.2f, like = %g.',
        #                 (iterCount/totIters)*100,
        #                 dbHatStar[2], dbHatStar[1], max(like, na.rm = TRUE) ))
      } else if (opt$progress == 2) {
        message(sprintf('Progress : %5.1f%%, b = %4.2f, d = %4.2f, like = %g.',
                        (iterCount/totIters)*100,
                        dbHatStar[2], dbHatStar[1], max(like, na.rm = TRUE) ))
      }

    }

  }


  #--------------------------------------------------------------------------------
  # STORE THE PARAMETER VALUES
  #--------------------------------------------------------------------------------

  # Output a list of values, including the previous vector \code{params}.
  likeGrid_params <- list(
    params = NA,
    dbHatStar = dbHatStar,
    muHatStar = NA,
    # max_like = max(like, na.rm = TRUE),
    max_like = max_like,
    local_max = NA,
    Grid2d = Grid2d,
    dGrid = dGrid,
    bGrid = bGrid,
    dGrid_orig = dGrid_orig,
    bGrid_orig = bGrid_orig,
    like = like
  )

  # Record local maxima, if required.
  if(opt$LocalMax) {
    likeGrid_params$local_max = local_max
  }


  # Set parameter values for output.
  params <- dbHatStar
  likeGrid_params$params <- params

  # Add level parameter corresponding to max likelihood.
  if(opt$levelParam) {

    muHatStar  <- mu[, indexB, indexD]
    likeGrid_params$muHatStar <- muHatStar


    # Concatenate parameter values.
    params <- matrix(c(params, muHatStar),
                     nrow = 1, ncol = length(params) + length(muHatStar))


    likeGrid_params$params <- params

  }

  # Append additional parameters and set class of output.
  likeGrid_params$k <- k
  likeGrid_params$r <- r
  likeGrid_params$opt <- opt
  class(likeGrid_params) <- 'FCVAR_grid'


  #--------------------------------------------------------------------------------
  # PLOT THE LIKELIHOOD
  #--------------------------------------------------------------------------------


  if (opt$plotLike) {

    graphics::plot(x = likeGrid_params)

  }


  return(likeGrid_params)
}


#' Plot the Likelihood Function for the FCVAR Model
#'
#' \code{plot.FCVAR_grid} plots the likelihood function from \code{FCVARlikeGrid}.
#' \code{FCVARlikeGrid} performs a grid-search optimization
#' by calculating the likelihood function
#' on a grid of candidate parameter values.
#' This function evaluates the likelihood over a grid of values
#' 	for \code{c(d,b)} (or \code{phi}, when  there are constraints on \code{c(d,b)}).
#' 	It can be used when parameter estimates are sensitive to
#' 	starting values to give an approximation of the global max which can
#' 	then be used as the starting value in the numerical optimization in
#' 	\code{FCVARestn}.
#'
#' @param x An S3 object of type \code{FCVAR_grid} output from \code{FCVARlikeGrid}.
#' @param y An argument for generic method \code{plot} that is not used in \code{plot.FCVAR_grid}.
#' @param ... Arguments to be passed to methods, such as graphical parameters
#' for the generic plot function.
#' @return NULL
#' @examples
#' \donttest{
#' opt <- FCVARoptions()
#' opt$dbStep1D     <- 0.1 # Coarser grid for plotting example.
#' opt$dbStep2D     <- 0.2 # Coarser grid for plotting example.
#' opt$gridSearch   <- 0 # Disable grid search in optimization.
#' opt$dbMin        <- c(0.01, 0.01) # Set lower bound for d,b.
#' opt$dbMax        <- c(2.00, 2.00) # Set upper bound for d,b.
#' opt$constrained  <- 0 # Impose restriction dbMax >= d >= b >= dbMin ? 1 <- yes, 0 <- no.
#' x <- votingJNP2014[, c("lib", "ir_can", "un_can")]
#' opt$progress <- 2 # Show progress report on each value of b.
#' likeGrid_params <- FCVARlikeGrid(x, k = 2, r = 1, opt)
#' graphics::plot(likeGrid_params)
#' }
#' @note Calls \code{graphics::persp} when \code{x$Grid2d == TRUE} and
#' calls \code{graphics::plot} when \code{x$Grid2d == FALSE}.
#' @family FCVAR auxiliary functions
#' @seealso \code{FCVARoptions} to set default estimation options.
#' \code{plot.FCVAR_grid} plots the likelihood function from \code{FCVARlikeGrid}.
#' @export
#'
plot.FCVAR_grid <- function(x, y = NULL, ...) {

  # Extract parameters.
  Grid2d <- x$Grid2d
  like <- x$like
  dGrid <- x$dGrid
  bGrid <- x$bGrid
  dGrid_orig <- x$dGrid_orig
  bGrid_orig <- x$bGrid_orig
  opt <- x$opt

  # Additional parameters.
  dots <- list(...)


  # if (!is.null(main)) {
  if ('main' %in% names(dots)) {
    main <- dots$main
  } else {
      main <- c('Log-likelihood Function ',
                sprintf('Rank: %d, Lags: %d', x$r, x$k))
  }


  # Plot likelihood depending on dimension of search.
  if(Grid2d) {
    # 2-dimensional plot.

    colors <- grDevices::rainbow(100)
    like2D <- t(like)
    like.facet.center <- (like2D[-1, -1] + like2D[-1, -ncol(like2D)] + like2D[-nrow(like2D), -1] + like2D[-nrow(like2D), -ncol(like2D)])/4
    # Range of the facet center on a 100-scale (number of colors)
    like.facet.range <- cut(like.facet.center, 100)

    # Reduce the size of margins.
    # First, store user's existing par() settings and restore on exit.
    oldpar <- graphics::par(no.readonly = TRUE)
    on.exit(graphics::par(oldpar))

    # Defaults are as follows:
    # graphics::par()$mar
    # 5.1 4.1 4.1 2.1
    # bottom, left, top and right margins respectively.
    graphics::par(mar = c(1.1, 1.1, 2.1, 1.1))

    # Remember to reset after:
    # Not with defaults:
    # graphics::par(mar = c(5.1, 4.1, 4.1, 2.1))
    # But with whatever user settings were in place:
    # graphics::par(oldpar)

    graphics::persp(dGrid_orig, bGrid_orig,
          like2D,
          xlab = 'd',
          ylab = 'b',
          zlab = '',
          main = main,
          phi = 45, theta = -55,
          r = 1.5,
          d = 4.0,
          xaxs = "i",
          col = colors[like.facet.range]
    )

    # Reset the size of margins.
    # graphics::par(mar = c(5.1, 4.1, 4.1, 2.1))
    graphics::par(oldpar)

  } else {
    # 1-dimensional plot.
    if ((opt$R_psi[1] == 1) & (opt$R_psi[2] == -1) & (opt$r_psi[1] == 0)) {
      like_x_label <- 'd = b'
      plot_grid <- bGrid_orig
    } else {
      like_x_label <- 'phi'
      plot_grid <- bGrid
    }

    graphics::plot(plot_grid, like,
         main = main,
         ylab = 'log-likelihood',
         xlab = like_x_label,
         type = 'l',
         col = 'blue',
         lwd = 3)

  }


}


# Find local optima
#
# \code{find_local_max} finds the row and column numbers corresponding
# to the local maxima of a 1- or 2-dimensional array.
# Note that Roxygen comments are excluded to keep this function internal.
# However, the contents of Roxygen comments are shown below for those who read the scripts.
#
# @param x A 1- or 2-dimensional numerical array of the candidate values
# of the objective function.
# @return A list object with three elements:
# \describe{
#   \item{\code{row}}{An integer row number of the local maxima. }
#   \item{\code{col}}{An integer column number of the local maxima. }
#   \item{\code{value}}{The numeric values of the local maxima. }
# }
#
#
find_local_max <- function(x) {

  # Determine size of input.
  nrows <- dim(x)[1]
  ncols <- dim(x)[2]

  # Remove corner cases.
  if (nrows == 1 & ncols == 1) {
    return(list(row = 1, col = 1, value = x))
  } else if (nrows == 1) {
    x <- t(x)
    nrows <- ncols
    ncols <- 1
  }

  # Initialize matrix of local max indicators.
  loc_max_ind <- matrix(TRUE, nrow = nrows, ncol = ncols)

  # Proceed by ruling out locally dominated points.

  # Check above.
  x_check <- matrix( - Inf, nrow = nrows, ncol = ncols)
  x_check[ - 1, ] <- x[ - nrows, ]
  loc_max_ind <- loc_max_ind & (x > x_check)

  # Check below.
  x_check <- matrix( - Inf, nrow = nrows, ncol = ncols)
  x_check[ - nrows, ] <- x[ - 1, ]
  loc_max_ind <- loc_max_ind & (x > x_check)


  # Check sides and diagonals only if 2-dimensional matrix.
  if (ncols > 1) {

    # Check left.
    x_check <- matrix( - Inf, nrow = nrows, ncol = ncols)
    x_check[ , - 1] <- x[ , - ncols]
    loc_max_ind <- loc_max_ind & (x > x_check)

    # Check right.
    x_check <- matrix( - Inf, nrow = nrows, ncol = ncols)
    x_check[ , - ncols] <- x[ , - 1]
    loc_max_ind <- loc_max_ind & (x > x_check)


    # Check top left.
    x_check <- matrix( - Inf, nrow = nrows, ncol = ncols)
    x_check[ - 1, - 1] <- x[ - nrows, - ncols]
    loc_max_ind <- loc_max_ind & (x > x_check)

    # Check top right.
    x_check <- matrix( - Inf, nrow = nrows, ncol = ncols)
    x_check[ - 1, - ncols] <- x[ - nrows, - 1]
    loc_max_ind <- loc_max_ind & (x > x_check)


    # Check bottom left.
    x_check <- matrix( - Inf, nrow = nrows, ncol = ncols)
    x_check[ - nrows, - 1] <- x[ - 1, - ncols]
    loc_max_ind <- loc_max_ind & (x > x_check)

    # Check bottom right.
    x_check <- matrix( - Inf, nrow = nrows, ncol = ncols)
    x_check[ - nrows, - ncols] <- x[ - 1, - 1]
    loc_max_ind <- loc_max_ind & (x > x_check)

  }


  # Assemble output into a list.
  loc_max_out <- list(
    row = matrix(rep(seq(nrows), ncols),
                 nrow = nrows,
                 ncol = ncols)[loc_max_ind],
    col = matrix(rep(seq(ncols), nrows),
                 nrow = nrows,
                 ncol = ncols, byrow = TRUE)[loc_max_ind],
    value = x[loc_max_ind])

  return(loc_max_out)
}



# Likelihood Function for the Unconstrained FCVAR Model
#
# \code{FCVARlikeMu} calculates the likelihood for the unconstrained FCVAR model
# for a given set of parameter values. It is used by the \code{FCVARlikeGrid}
# function to numerically optimize over the level parameter for given values of
# the fractional parameters.
# Note that Roxygen comments are excluded to keep this function internal.
# However, the contents of Roxygen comments are shown below for those who read the scripts.
#
# @param y A matrix of variables to be included in the system.
# @param db The orders of fractional integration.
# @param mu The level parameter.
# @param k The number of lags in the system.
# @param r The cointegrating rank.
# @param opt An S3 object of class \code{FCVAR_opt} that stores the chosen estimation options,
# generated from \code{FCVARoptions()}.
# @return A number \code{like}, the log-likelihood evaluated at the
# specified parameter values.
# @examples
# opt <- FCVARoptions()
# x <- votingJNP2014[, c("lib", "ir_can", "un_can")]
# like <- FCVARlikeMu(mu = colMeans(x), y = x, db = c(1, 1), k = 2, r = 1, opt)
# @family FCVAR auxiliary functions
# @seealso \code{FCVARoptions} to set default estimation options.
# The \code{FCVARlikeGrid} calls this function to perform a grid search over the
# parameter values.
# @export
#
FCVARlikeMu <- function(mu, y, db, k, r, opt) {


  t <- nrow(y)
  x <- y - matrix(1, nrow = t, ncol = 1) %*% mu


  # Obtain concentrated parameter estimates.
  estimates <- GetParams(x, k, r, db, opt)

  # Calculate value of likelihood function.
  cap_T <- t - opt$N
  p <- ncol(x)
  like <- - cap_T*p/2*( log(2*pi) + 1)  - cap_T/2*log(det(estimates$OmegaHat))

  return(like)
}


# Likelihood Function for the Constrained FCVAR Model
#
# \code{FCVARlike} calculates the likelihood for the constrained FCVAR model
# for a given set of parameter values. This function adjusts the variables with the level parameter,
# if present, and returns the log-likelihood given \code{d} and \code{b}.
# Note that Roxygen comments are excluded to keep this function internal.
# However, the contents of Roxygen comments are shown below for those who read the scripts.
#
# @param params A vector of parameters \code{d} and \code{b}
# (and \code{mu} if option selected).
# @param x A matrix of variables to be included in the system.
# @param k The number of lags in the system.
# @param r The cointegrating rank.
# @param opt An S3 object of class \code{FCVAR_opt} that stores the chosen estimation options,
# generated from \code{FCVARoptions()}.
# @return A number \code{like}, the log-likelihood evaluated at the
# specified parameter values.
# @examples
# opt <- FCVARoptions()
# opt$gridSearch   <- 0 # Disable grid search in optimization.
# opt$dbMin        <- c(0.01, 0.01) # Set lower bound for d,b.
# opt$dbMax        <- c(2.00, 2.00) # Set upper bound for d,b.
# opt$constrained  <- 0 # Impose restriction dbMax >= d >= b >= dbMin ? 1 <- yes, 0 <- no.
# x <- votingJNP2014[, c("lib", "ir_can", "un_can")]
# results <- FCVARestn(x, k = 2, r = 1, opt)
# FCVARlike(c(results$coeffs$db, results$coeffs$muHat), x, k = 2, r = 1, opt)
# @family FCVAR auxiliary functions
# @seealso \code{FCVARoptions} to set default estimation options.
# @export
#
FCVARlike <- function(params, x, k, r, opt) {

  # global estimatesTEMP

  cap_T <- nrow(x)
  p <- ncol(x)


  # If there is one linear restriction imposed, optimization is over phi,
  #  so translate from phi to (d,b), otherwise, take parameters as
  #  given.
  if(!is.null((opt$R_psi)) && nrow(opt$R_psi) == 1) {
    H <- pracma::null(opt$R_psi)
    h <- t(opt$R_psi) %*% solve(opt$R_psi %*% t(opt$R_psi)) %*% opt$r_psi
    db <- H*params[1] + h
    d <- db[1]
    b <- db[2]
    if (opt$levelParam) {
      mu <- params[2:length(params)]
    }

  } else {
    d <- params[1]
    b <- params[2]
    if (opt$levelParam) {
      mu <- params[3:length(params)]
    }

  }


  # Modify the data by a mu level.
  if (opt$levelParam) {

    y <- x - matrix(1, nrow = cap_T, ncol = p) %*% diag(mu)
  } else {
    y <- x
  }


  # If d=b is not imposed via opt$restrictDB, but d>=b is required in
  #	opt$constrained, then impose the inequality here.
  if(opt$constrained) {
    b <- min(b, d)
  }


  # Obtain concentrated parameter estimates.
  estimates <- GetParams(y, k, r, c(d, b), opt)



  # Storing the estimates in a global structure here allows us to skip a
  #   call to GetParams after optimization to recover the coefficient
  #   estimates
  estimatesTEMP <- estimates
  # If level parameter is present, they are the last p parameters in the
  #   params vector
  if (opt$levelParam) {
    estimatesTEMP$muHat <- mu
  } else {
    estimatesTEMP$muHat <- NULL
  }




  # Calculate value of likelihood function.
  # Base this on raw inputs:
  cap_T <- nrow(x) - opt$N
  p <- ncol(x)
  like <- - cap_T*p/2*( log(2*pi) + 1)  - cap_T/2*log(det(estimates$OmegaHat))


  return(like)
}


# Obtain concentrated parameter estimates for the Constrained FCVAR Model
#
# \code{GetEstimates} calculates the concentrated parameter estimates for
# the constrained FCVAR model for a given set of parameter values.
# It is used after estimating the parameters that are numerically optimized
# to obtain the corresponding concentrated parameters.
# Like \code{FCVARlike} , this function adjusts the variables with the level parameter,
# if present, and returns the log-likelihood given \code{d} and \code{b}.
# Note that Roxygen comments are excluded to keep this function internal.
# However, the contents of Roxygen comments are shown below for those who read the scripts.
#
# @param params A vector of parameters \code{d} and \code{b}
# (and \code{mu} if option selected).
# @param x A matrix of variables to be included in the system.
# @param k The number of lags in the system.
# @param r The cointegrating rank.
# @param opt An S3 object of class \code{FCVAR_opt} that stores the chosen estimation options,
# generated from \code{FCVARoptions()}.
# @return A list object \code{estimates} containing the estimates,
# including the following parameters:
# \describe{
#   \item{\code{db}}{The orders of fractional integration (taken directly from the input).}
#   \item{\code{alphaHat}}{A \eqn{p x r} matrix of adjustment parameters.}
#   \item{\code{betaHat}}{A \eqn{p x r} matrix of cointegrating vectors.
#       The \eqn{r x 1} vector \eqn{\beta x_t} is the stationary cointegration relations.}
#   \item{\code{rhoHat}}{A \eqn{p x 1} vector of restricted constatnts.}
#   \item{\code{piHat}}{A \eqn{p x p} matrix \eqn{\Pi = \alpha \beta'} of long-run levels.}
#   \item{\code{OmegaHat}}{A \eqn{p x p} covariance matrix of the error terms.}
#   \item{\code{GammaHat}}{A ( \eqn{p x kp} matrix \code{cbind(GammaHat1,...,GammaHatk)})
#   of autoregressive coefficients. }
#   \item{\code{muHat}}{A vector of the optimal \code{mu} if level parameter is selected. }
# }
# @examples
# opt <- FCVARoptions()
# opt$gridSearch   <- 0 # Disable grid search in optimization.
# opt$dbMin        <- c(0.01, 0.01) # Set lower bound for d,b.
# opt$dbMax        <- c(2.00, 2.00) # Set upper bound for d,b.
# opt$constrained  <- 0 # Impose restriction dbMax >= d >= b >= dbMin ? 1 <- yes, 0 <- no.
# x <- votingJNP2014[, c("lib", "ir_can", "un_can")]
# results <- FCVARestn(x, k = 2, r = 1, opt)
# GetEstimates(c(results$coeffs$db, results$coeffs$muHat), x, k = 2, r = 1, opt)
# @family FCVAR auxiliary functions
# @seealso \code{FCVARoptions} to set default estimation options.
# \code{FCVARlike} performs the same calculations to obtain the value
# of the likelihood function.
# @export
#
GetEstimates <- function(params, x, k, r, opt) {

  cap_T <- nrow(x)
  p <- ncol(x)


  # If there is one linear restriction imposed, optimization is over phi,
  #  so translate from phi to (d,b), otherwise, take parameters as
  #  given.
  if(!is.null((opt$R_psi)) && nrow(opt$R_psi) == 1) {
    H <- pracma::null(opt$R_psi)
    h <- t(opt$R_psi) %*% solve(opt$R_psi %*% t(opt$R_psi)) %*% opt$r_psi
    db <- H*params[1] + h
    d <- db[1]
    b <- db[2]
    if (opt$levelParam) {
      mu <- params[2:length(params)]
    }

  } else {
    d <- params[1]
    b <- params[2]
    if (opt$levelParam) {
      mu <- params[3:length(params)]
    }

  }


  # Modify the data by a mu level.
  if (opt$levelParam) {

    y <- x - matrix(1, nrow = cap_T, ncol = p) %*% diag(mu)
  } else {
    y <- x
  }


  # If d=b is not imposed via opt$restrictDB, but d>=b is required in
  #	opt$constrained, then impose the inequality here.
  if(opt$constrained) {
    b <- min(b, d)
  }


  # Obtain concentrated parameter estimates$
  estimates <- GetParams(y, k, r, c(d, b), opt)

  # If level parameter is present, they are the last p parameters in the
  #   params vector
  if (opt$levelParam) {
    estimates$muHat <- mu
  } else {
    estimates$muHat <- NULL
  }

  return(estimates)

}


# Likelihood Function for the FCVAR Model
#
# \code{FCVARlikeFull} calculates the likelihood for the constrained FCVAR model
# for a given set of parameter values.
# This function takes the full set of coefficients from a call
# to \code{FCVARestn} and returns the log-likelihood given \code{d} and \code{b}.
#
# @param x A matrix of variables to be included in the system.
# @param k The number of lags in the system.
# @param r The cointegrating rank.
# @param coeffs A list of coefficients of the FCVAR model.
# It is an element of the list \code{results} returned by \code{FCVARestn},
# without the parameters \code{betaHat} and \code{rhoHat}.
# @param betaHat A \eqn{p x r} matrix of cointegrating vectors.
# The \eqn{r x 1} vector \eqn{\beta x_t} is the stationary cointegration relations.
# @param rhoHat A \eqn{p x 1} vector of restricted constatnts.
# @param opt An S3 object of class \code{FCVAR_opt} that stores the chosen estimation options,
# generated from \code{FCVARoptions()}.
# @return A number \code{like}, the log-likelihood evaluated at the
# specified parameter values.
# @examples
# opt <- FCVARoptions()
# opt$gridSearch   <- 0 # Disable grid search in optimization.
# opt$dbMin        <- c(0.01, 0.01) # Set lower bound for d,b.
# opt$dbMax        <- c(2.00, 2.00) # Set upper bound for d,b.
# opt$constrained  <- 0 # Impose restriction dbMax >= d >= b >= dbMin ? 1 <- yes, 0 <- no.
# x <- votingJNP2014[, c("lib", "ir_can", "un_can")]
# results <- FCVARestn(x, k = 2, r = 1, opt)
# FCVARlikeFull(x, k = 2, r = 1, coeffs = results$coeffs,
#               beta = results$coeffs$betaHat, rho = results$coeffs$rhoHat, opt)
# @family FCVAR auxiliary functions
# @seealso \code{FCVARoptions} to set default estimation options.
# \code{FCVARestn} for the estimation of coefficients in \code{coeffs}.
# @export
#
FCVARlikeFull <- function(x, k, r, coeffs, betaHat, rhoHat, opt) {


  # Add betaHat and rhoHat to the coefficients to get residuals because they
  #   are not used in the Hessian calculation and are missing from the
  #   structure coeffs
  coeffs$betaHat <- betaHat
  coeffs$rhoHat <- rhoHat

  # Obtain residuals.
  epsilon <- GetResiduals(x, k, r, coeffs, opt)


  # Calculate value of likelihood function.
  cap_T <- nrow(x) - opt$N
  p <- ncol(x)
  OmegaHat <- t(epsilon) %*% epsilon/cap_T
  like <- - cap_T*p/2*( log(2*pi) + 1)  - cap_T/2*log(det(OmegaHat))

  return(like)
}


# Transform Data for Regression
#
# \code{TransformData} transforms the raw data by fractional differencing.
# The output is in the format required for regression and
# reduced rank regression.
# Note that Roxygen comments are excluded to keep this function internal.
# However, the contents of Roxygen comments are shown below for those who read the scripts.
#
# @param x A matrix of variables to be included in the system.
# @param k The number of lags in the system.
# @param db The orders of fractional integration.
# @param opt An S3 object of class \code{FCVAR_opt} that stores the chosen estimation options,
# generated from \code{FCVARoptions()}.
# @return A list object \code{Z_array} containing the transformed data,
# including the following parameters:
# \describe{
#   \item{\code{Z0}}{A matrix of data calculated by fractionally differencing
#   \code{x} at differencing order \code{d}.}
#   \item{\code{Z1}}{The matrix \code{x} (augmented with a vector of ones
#   if the model includes a restricted constant term), which is then lagged and stacked
#   and fractionally differenced (with order \eqn{d-b}).}
#   \item{\code{Z2}}{The matrix \code{x} lagged and stacked
#   and fractionally differenced (with order \code{d}).}
#   \item{\code{Z3}}{A column of ones if model includes an
#   unrestricted constant term, otherwise \code{NULL}.}
# }
# @examples
# opt <- FCVARoptions()
# opt$gridSearch   <- 0 # Disable grid search in optimization.
# opt$dbMin        <- c(0.01, 0.01) # Set lower bound for d,b.
# opt$dbMax        <- c(2.00, 2.00) # Set upper bound for d,b.
# opt$constrained  <- 0 # Impose restriction dbMax >= d >= b >= dbMin ? 1 <- yes, 0 <- no.
# x <- votingJNP2014[, c("lib", "ir_can", "un_can")]
# results <- FCVARestn(x, k = 2, r = 1, opt)
# Z_array <- TransformData(x, k = 2, db = results$coeffs$db, opt)
# @family FCVAR auxiliary functions
# @seealso \code{FCVARoptions} to set default estimation options.
# \code{FCVARestn} calls \code{GetParams}, which calls \code{TransformData}
# to estimate the FCVAR model.
# \code{TransformData} in turn calls \code{FracDiff} and \code{Lbk}
# to perform the transformation.
# @references Johansen, S. and M. \enc{Ø}{O}. Nielsen (2012).
# "Likelihood inference for a fractionally cointegrated
# vector autoregressive model," Econometrica 80, 2667-2732.
# @references Johansen, S. (1995). "Likelihood-Based Inference
# in Cointegrated Vector Autoregressive Models,"
# New York: Oxford University Press.
# @export
#
TransformData <- function(x, k, db, opt) {

  # Number of initial values and sample size.
  N <- opt$N
  cap_T <- nrow(x) - N

  # Extract parameters from input.
  d <- db[1]
  b <- db[2]

  # Transform data as required.
  Z0 <- FracDiff(x, d)

  Z1 <- x
  # Add a column with ones if model includes a restricted constant term.
  if(opt$rConstant) {
    Z1 <- cbind(x, matrix(1, nrow = N + cap_T, ncol = 1))
  }

  Z1 <- FracDiff(  Lbk( Z1 , b, 1)  ,  d - b )

  Z2 <- FracDiff(  Lbk( x , b, k)  , d)

  # Z3 contains the unrestricted deterministics
  Z3 <- NULL
  # # Add a column with ones if model includes a unrestricted constant term.
  if(opt$unrConstant) {
    Z3 <- matrix(1, nrow = cap_T, ncol = 1)
  }


  # Trim off initial values.
  Z0 <- Z0[(N+1):nrow(Z0), ]
  Z1 <- Z1[(N+1):nrow(Z1), ]
  if(k > 0) {
    Z2 <- Z2[(N+1):nrow(Z2), ]
  }




  # Return all arrays together as a list.
  Z_array <- list(
    Z0 = Z0,
    Z1 = Z1,
    Z2 = Z2,
    Z3 = Z3
  )

  return(Z_array)
}


# Calculate Residuals for the FCVAR Model
#
# \code{GetResiduals} calculates residuals for the FCVAR model
# from given parameter values.
# Note that Roxygen comments are excluded to keep this function internal.
# However, the contents of Roxygen comments are shown below for those who read the scripts.
#
# @param x A matrix of variables to be included in the system.
# @param k The number of lags in the system.
# @param r The cointegrating rank.
# @param coeffs A list of coefficients of the FCVAR model.
# It is an element of the list \code{results} returned by \code{FCVARestn}.
# @param opt An S3 object of class \code{FCVAR_opt} that stores the chosen estimation options,
# generated from \code{FCVARoptions()}.
# @return A matrix \code{epsilon} of residuals from FCVAR model estimation
# calculated with the parameter estimates specified in \code{coeffs}.
# @examples
# opt <- FCVARoptions()
# opt$gridSearch   <- 0 # Disable grid search in optimization.
# opt$dbMin        <- c(0.01, 0.01) # Set lower bound for d,b.
# opt$dbMax        <- c(2.00, 2.00) # Set upper bound for d,b.
# opt$constrained  <- 0 # Impose restriction dbMax >= d >= b >= dbMin ? 1 <- yes, 0 <- no.
# x <- votingJNP2014[, c("lib", "ir_can", "un_can")]
# results <- FCVARestn(x, k = 2, r = 1, opt)
# epsilon <- GetResiduals(x, k = 2, r = 1, coeffs = results$coeffs, opt)
# @family FCVAR auxiliary functions
# @seealso \code{FCVARoptions} to set default estimation options.
# \code{FCVARestn} to estimate the FCVAR model.
# @references Johansen, S. and M. \enc{Ø}{O}. Nielsen (2012).
# "Likelihood inference for a fractionally cointegrated
# vector autoregressive model," Econometrica 80, 2667-2732.
# @export
#
GetResiduals <- function(x, k, r, coeffs, opt) {


  # If level parameter is included, the data must be shifted before
  #   calculating the residuals:
  if (opt$levelParam) {
    cap_T <- nrow(x)
    y <- x - matrix(1, nrow = cap_T, ncol = 1) %*% coeffs$muHat
  } else {
    y <- x
  }

  #--------------------------------------------------------------------------------
  # Transform data
  #--------------------------------------------------------------------------------
  Z_array <- TransformData(y, k, coeffs$db, opt)
  Z0 <- Z_array$Z0
  Z1 <- Z_array$Z1
  Z2 <- Z_array$Z2
  Z3 <- Z_array$Z3

  #--------------------------------------------------------------------------------
  # Calculate residuals
  #--------------------------------------------------------------------------------
  epsilon <- Z0

  if (r > 0) {
    epsilon <- epsilon -
      Z1 %*% rbind(coeffs$betaHat, coeffs$rhoHat) %*% t(coeffs$alphaHat)
  }


  if (k > 0) {
    epsilon <- epsilon - Z2 %*% t(coeffs$GammaHat)
  }


  if (opt$unrConstant) {
    epsilon <- epsilon - Z3 %*% t(coeffs$xiHat)
  }



  return(epsilon)
}


# Calculate Lag Polynomial in the Fractional Lag Operator
#
# \code{Lbk} calculates a lag polynomial in the fractional lag operator.
# Note that Roxygen comments are excluded to keep this function internal.
# However, the contents of Roxygen comments are shown below for those who read the scripts.
#
# @param x A matrix of variables to be included in the system.
# @param b The order of fractional differencing.
# @param k The number of lags in the system.
# @return A matrix \code{Lbkx} of the form \eqn{[ Lb^1 x, Lb^2 x, ..., Lb^k x]}
# where \eqn{Lb = 1 - (1-L)^b}.
# The output matrix has the same number of rows as \code{x}
# but \code{k} times as many columns.
# @examples
# opt <- FCVARoptions()
# opt$gridSearch   <- 0 # Disable grid search in optimization.
# opt$dbMin        <- c(0.01, 0.01) # Set lower bound for d,b.
# opt$dbMax        <- c(2.00, 2.00) # Set upper bound for d,b.
# opt$constrained  <- 0 # Impose restriction dbMax >= d >= b >= dbMin ? 1 <- yes, 0 <- no.
# x <- votingJNP2014[, c("lib", "ir_can", "un_can")]
# results <- FCVARestn(x, k = 2, r = 1, opt)
# Lbkx <- Lbk(x, b = results$coeffs$db[2], k = 2)
# @family FCVAR auxiliary functions
# @seealso \code{FCVARoptions} to set default estimation options.
# \code{FCVARestn} calls \code{GetParams}, which calls \code{TransformData}
# to estimate the FCVAR model.
# \code{TransformData} in turn calls \code{FracDiff} and \code{Lbk}
# to perform the transformation.
# @export
#
Lbk <- function(x, b, k) {

  p <- ncol(x)

  # Initialize output matrix.
  Lbkx <- NULL

  # For i <- 1, set first column of Lbkx <- Lb^1 x.
  if (k > 0) {
    bx <- FracDiff(x, b)
    Lbkx <- x - bx
  }


  if (k > 1) {
    for (i in 2:k ) {
      Lbkx <- cbind(Lbkx,
                    ( Lbkx[ , (p*(i-2)+1) : ncol(Lbkx)] -
                        FracDiff(Lbkx[ , (p*(i-2)+1) : ncol(Lbkx)], b) ))

    }

  }



  return(Lbkx)
}



#' Fast Fractional Differencing
#'
#' \code{FracDiff} is a fractional differencing procedure based on the
#' 	fast fractional difference algorithm of Jensen & Nielsen (2014).
#'
#' @param x A matrix of variables to be included in the system.
#' @param d The order of fractional differencing.
#' @return A vector or matrix \code{dx} equal to \eqn{(1-L)^d x}
#' of the same dimensions as x.
#' @examples
#' set.seed(42)
#' WN <- matrix(stats::rnorm(200), nrow = 100, ncol = 2)
#' \donttest{MVWNtest_stats <- MVWNtest(x = WN, maxlag = 10, printResults = 1)}
#' x <- FracDiff(x = WN, d = - 0.5)
#' \donttest{MVWNtest_stats <- MVWNtest(x = x, maxlag = 10, printResults = 1)}
#' WN_x_d <- FracDiff(x, d = 0.5)
#' \donttest{MVWNtest_stats <- MVWNtest(x = WN_x_d, maxlag = 10, printResults = 1)}
#' @family FCVAR auxiliary functions
#' @seealso \code{FCVARoptions} to set default estimation options.
#' \code{FCVARestn} calls \code{GetParams}, which calls \code{TransformData}
#' to estimate the FCVAR model.
#' \code{TransformData} in turn calls \code{FracDiff} and \code{Lbk}
#' to perform the transformation.
#' @references Jensen, A. N. and M. \enc{Ø}{O}. Nielsen (2014).
#' "A fast fractional difference algorithm,"
#' Journal of Time Series Analysis 35, 428-436.
#' @note This function differs from the \code{diffseries} function
#' in the \code{fracdiff} package, in that the \code{diffseries}
#' function demeans the series first.
#' In particular, the difference between the out put of the function calls
#' \code{FCVAR::FracDiff(x - mean(x), d = 0.5)}
#' and \code{fracdiff::diffseries(x, d = 0.5)} is numerically small.
#' @export
#'
FracDiff <- function(x, d) {


  if(is.null(x)) {
    dx <- NULL
  } else {


    p <- ncol(x)
    cap_T <- nrow(x)
    np2 <- stats::nextn(2*cap_T - 1, 2)

    k <- 1:(cap_T-1)
    b <- c(1, cumprod((k - d - 1)/k))


    dx <- matrix(0, nrow = cap_T, ncol = p)

    for (i in 1:p) {


      dxi <- stats::fft(stats::fft(c(b, rep(0, np2 - cap_T))) *
                          stats::fft(c(x[, i], rep(0, np2 - cap_T))), inverse = cap_T) / np2

      dx[, i] <- Re(dxi[1:cap_T])

    }



  }

  return(dx)
}


# Count the Number of Free Parameters
#
# \code{GetFreeParams} counts the number of free parameters based on
# 	the number of coefficients to estimate minus the total number of
# 	restrictions. When both \code{alpha} and \code{beta} are restricted,
# 	the rank condition is used to count the free parameters in those two variables.
# Note that Roxygen comments are excluded to keep this function internal.
# However, the contents of Roxygen comments are shown below for those who read the scripts.
#
# @param k The number of lags in the system.
# @param r The cointegrating rank.
# @param p The number of variables in the system.
# @param opt An S3 object of class \code{FCVAR_opt} that stores the chosen estimation options,
# generated from \code{FCVARoptions()}.
# @param rankJ The rank of a conditioning matrix, as described in
# Boswijk & Doornik (2004, p.447), which is only used if there are
# restrictions imposed on \code{alpha} or \code{beta}, otherwise \code{NULL}.
# @return The number of free parameters \code{fp}.
# @examples
# opt <- FCVARoptions()
# opt$gridSearch   <- 0 # Disable grid search in optimization.
# opt$dbMin        <- c(0.01, 0.01) # Set lower bound for d,b.
# opt$dbMax        <- c(2.00, 2.00) # Set upper bound for d,b.
# opt$constrained  <- 0 # Impose restriction dbMax >= d >= b >= dbMin ? 1 <- yes, 0 <- no.
# newOpt <- FCVARoptionUpdates(opt, p = 3, r = 1) # Need to update restriction matrices.
# GetFreeParams(k = 2, r = 1, p = 3, opt = newOpt, rankJ = NULL)
#
# opt <- FCVARoptions()
# opt$dbMin        <- c(0.01, 0.01) # Set lower bound for d,b.
# opt$dbMax        <- c(2.00, 2.00) # Set upper bound for d,b.
# opt$constrained  <- 0 # Impose restriction dbMax >= d >= b >= dbMin ? 1 <- yes, 0 <- no.
# opt$R_Beta <- matrix(c(1, 0, 0), nrow = 1, ncol = 3)
# newOpt <- FCVARoptionUpdates(opt, p = 3, r = 1) # Need to update restriction matrices.
# GetFreeParams(k = 2, r = 1, p = 3, opt = newOpt, rankJ = 4)
#
# opt <- FCVARoptions()
# opt$dbMin        <- c(0.01, 0.01) # Set lower bound for d,b.
# opt$dbMax        <- c(2.00, 2.00) # Set upper bound for d,b.
# opt$constrained  <- 0 # Impose restriction dbMax >= d >= b >= dbMin ? 1 <- yes, 0 <- no.
# opt$R_Alpha <- matrix(c(0, 1, 0), nrow = 1, ncol = 3)
# newOpt <- FCVARoptionUpdates(opt, p = 3, r = 1) # Need to update restriction matrices.
# GetFreeParams(k = 2, r = 1, p = 3, opt = newOpt, rankJ = 4)
# @family FCVAR auxiliary functions
# @seealso \code{FCVARoptions} to set default estimation options.
# \code{FCVARestn}, \code{HypoTest} and \code{LagSelect} to estimate the FCVAR model
# and use this in the calculation of the degrees of freedom
# for a variety of statistics.
# @references Boswijk, H. P. and J. A. Doornik (2004).
# "Identifying, estimating and testing restricted cointegrated systems:
# An overview," Statistica Neerlandica 58, 440-465.
# @export
#
GetFreeParams <- function(k, r, p, opt, rankJ) {

  # ---- First count the number of parameters -------- %
  fDB <- 2 # updated by rDB below for the model d=b (1 fractional parameter)
  fpA <- p*r # alpha
  fpB <- p*r # beta
  fpG <- p*p*k # Gamma
  fpM <- p*opt$levelParam # mu
  fpRrh <- r*opt$rConstant # restricted constant
  fpUrh <- p*opt$unrConstant # unrestricted constant

  numParams <- fDB + fpA + fpB + fpG + fpM + fpRrh + fpUrh

  # ---- Next count the number of restrictions -------- %
  # rDB <- nrow(opt$R_psi) # Returns numeric(0) when d==b not imposed.
  if(is.null(opt$R_psi)) {
    rDB <- 0
  } else {
    rDB <- nrow(opt$R_psi)
  }


  if(is.null(opt$R_Beta) & is.null(opt$R_Alpha)) {
    # If Alpha or Beta are unrestricted, only an identification restriction
    #   is imposed
    rB <- r*r # identification restrictions (eye(r))
    numRest <- rDB + rB
    # ------ Free parameters is the difference between the two ---- %
    fp <- numParams - numRest
  }
  else {
    # If there are restrictions imposed on beta or alpha, we can use the rank
    #   condition calculated in the main estimation function to give us the
    #   number of free parameters in alpha, beta, and the restricted constant
    #   if it is being estimated:
    fpAlphaBetaRhoR <- rankJ
    fp <- fpAlphaBetaRhoR + (fDB + fpG + fpM + fpUrh) - rDB
  }


  return(fp)
}


# Calculate the Hessian Matrix
#
# \code{FCVARhess} calculates the Hessian matrix of the
# 	log-likelihood by taking numerical derivatives.
# 	It is used to calculate the standard errors of parameter estimates.
# Note that Roxygen comments are excluded to keep this function internal.
# However, the contents of Roxygen comments are shown below for those who read the scripts.
#
# @param x A matrix of variables to be included in the system.
# @param k The number of lags in the system.
# @param r The cointegrating rank.
# @param coeffs A list of coefficients of the FCVAR model.
# It is an element of the list \code{results} returned by \code{FCVARestn}.
# @param opt An S3 object of class \code{FCVAR_opt} that stores the chosen estimation options,
# generated from \code{FCVARoptions()}.
# @return The \code{hessian} matrix  of second derivatives of the FCVAR
# log-likelihood function, calculated with the parameter estimates
# specified in \code{coeffs}.
# @examples
# opt <- FCVARoptions()
# opt$gridSearch   <- 0 # Disable grid search in optimization.
# opt$dbMin        <- c(0.01, 0.01) # Set lower bound for d,b.
# opt$dbMax        <- c(2.00, 2.00) # Set upper bound for d,b.
# opt$constrained  <- 0 # Impose restriction dbMax >= d >= b >= dbMin ? 1 <- yes, 0 <- no.
# x <- votingJNP2014[, c("lib", "ir_can", "un_can")]
# results <- FCVARestn(x, k = 2, r = 1, opt)
# hessian <- FCVARhess(x, k = 2, r = 1, coeffs = results$coeffs, opt)
# @family FCVAR auxiliary functions
# @seealso \code{FCVARoptions} to set default estimation options.
# \code{FCVARestn} to estimate the FCVAR model and calculate
# standard errors of the estimates.
# @export
#
FCVARhess <- function(x, k, r, coeffs, opt) {


  # Set dimensions of matrices.
  p <- ncol(x)

  # rhoHat and betaHat are not used in the Hessian calculation.
  rho <- coeffs$rhoHat
  beta <- coeffs$betaHat

  # Specify delta (increment for numerical derivatives).
  # delta <- 10^(-1)
  # delta <- 10^(-2)
  # delta <- 10^(-3)
  # delta <- 10^(-4) # Default.
  # delta <- 10^(-5)
  # delta <- 10^(-6)
  # delta <- 10^(-8)
  # Added as a parameter in FCVAR_options.
  delta <- opt$hess_delta

  # Convert the parameters to vector form.
  phi0 <- SEmat2vecU(coeffs, k, r, p, opt)

  # Calculate vector of increments.
  deltaPhi <- delta*matrix(1, nrow = nrow(phi0), ncol = ncol(phi0))

  nPhi <- length(phi0)

  # Initialize the Hessian matrix.
  hessian <- matrix(0, nrow = nPhi, ncol = nPhi)

  # The loops below evaluate the likelihood at second order incremental
  #   shifts in the parameters.


  for (i in 1:nPhi) {


    for (j in 1:i) {


      # positive shift in both parameters.
      phi1 <- phi0 + deltaPhi*( (1:nPhi) == i ) + deltaPhi*( (1:nPhi) == j )
      coeffsAdj <- SEvec2matU( phi1, k, r, p, opt )
      # calculate likelihood
      like1 <- FCVARlikeFull(x, k, r, coeffsAdj, beta, rho, opt )

      # negative shift in first parameter, positive shift in second
      # parameter. If same parameter, no shift.
      phi2 <- phi0 - deltaPhi*( (1:nPhi) == i ) + deltaPhi*( (1:nPhi) == j )
      coeffsAdj <- SEvec2matU( phi2, k, r, p, opt )
      # calculate likelihood
      like2 <- FCVARlikeFull(x, k, r, coeffsAdj, beta, rho, opt)

      # positive shift in first parameter, negative shift in second
      # parameter. If same parameter, no shift.
      phi3 <- phi0 + deltaPhi*( (1:nPhi) == i ) - deltaPhi*( (1:nPhi) == j )
      coeffsAdj <- SEvec2matU( phi3, k, r, p, opt )
      # calculate likelihood
      like3 <- FCVARlikeFull(x, k, r, coeffsAdj, beta, rho, opt)

      # negative shift in both parameters.
      phi4 <- phi0 - deltaPhi*( (1:nPhi) == i ) - deltaPhi*( (1:nPhi) == j )
      coeffsAdj <- SEvec2matU( phi4, k, r, p, opt )
      # calculate likelihood
      like4 <- FCVARlikeFull(x, k, r, coeffsAdj, beta, rho, opt)

      # The numerical approximation to the second derivative.
      hessian[i,j] <- ( like1 - like2 - like3 + like4 )/4/deltaPhi[i]/deltaPhi[j]

      # Hessian is symmetric.
      hessian[j,i] <- hessian[i,j]


    }

  }

  return(hessian)
}


# Collect Parameters into a Vector
#
# \code{SEmat2vecU} transforms the model parameters in matrix
# 	form into a vector.
# Note that Roxygen comments are excluded to keep this function internal.
# However, the contents of Roxygen comments are shown below for those who read the scripts.
#
# @param coeffs A list of coefficients of the FCVAR model.
# It is an element of the list \code{results} returned by \code{FCVARestn}.
# @param k The number of lags in the system.
# @param r The cointegrating rank.
# @param p The number of variables in the system.
# @param opt An S3 object of class \code{FCVAR_opt} that stores the chosen estimation options,
# generated from \code{FCVARoptions()}.
# @return A vector \code{param} of parameters in the FCVAR model.
# @examples
# opt <- FCVARoptions()
# opt$gridSearch   <- 0 # Disable grid search in optimization.
# opt$dbMin        <- c(0.01, 0.01) # Set lower bound for d,b.
# opt$dbMax        <- c(2.00, 2.00) # Set upper bound for d,b.
# opt$constrained  <- 0 # Impose restriction dbMax >= d >= b >= dbMin ? 1 <- yes, 0 <- no.
# x <- votingJNP2014[, c("lib", "ir_can", "un_can")]
# results <- FCVARestn(x, k = 2, r = 1, opt)
# params <- SEmat2vecU(coeffs = results$coeffs, k = 2, r = 1, p = 3, opt)
# coeffs <- SEvec2matU(param = params, k = 2, r = 1, p = 3, opt )
#
# params <- matrix(seq(25))
# coeffs <- SEvec2matU(param = params, k = 2, r = 1, p = 3, opt )
# params <- SEmat2vecU(coeffs = coeffs, k = 2, r = 1, p = 3, opt)
# @family FCVAR auxiliary functions
# @seealso \code{FCVARoptions} to set default estimation options.
# \code{FCVARestn} to estimate the FCVAR model and calculate
# standard errors of the estimates.
# \code{SEmat2vecU} is called by \code{FCVARhess} to sort the parameters
# into a vector to calculate the Hessian matrix.
# \code{SEvec2matU} is a near inverse of \code{SEmat2vecU},
# in the sense that \code{SEvec2matU} obtains only a
# subset of the parameters in \code{results$coeffs}.
#
SEmat2vecU <- function(coeffs, k, r, p , opt) {

  # If restriction d=b is imposed, only adjust d.
  if (opt$restrictDB) {
    param <- coeffs$db[1]
  }
  else {
    param <- coeffs$db
  }


  # Append the level parameter MuHat.
  if (opt$levelParam) {
    param <- cbind(matrix(param, nrow = 1, ncol = length(param)),
                   matrix(coeffs$muHat, nrow = 1, ncol = p))
  }


  # Unrestricted constant.
  if (opt$unrConstant) {
    param <- cbind(matrix(param, nrow = 1, ncol = length(param)),
                   matrix(coeffs$xiHat, nrow = 1, ncol = p))
  }


  # alphaHat
  if (r > 0) {
    param <- cbind(matrix(param, nrow = 1, ncol = length(param)),
                   matrix( coeffs$alphaHat, nrow = 1, ncol = p*r ))
  }


  # GammaHat
  if (k > 0) {
    param <- cbind(matrix(param, nrow = 1, ncol = length(param)),
                   matrix( coeffs$GammaHat, nrow = 1, ncol = p*p*k ))
  }


  return(param)
}


# Extract Parameters from a Vector
#
# \code{SEvec2matU} transforms the vectorized model parameters
# 	into matrices.
# Note that Roxygen comments are excluded to keep this function internal.
# However, the contents of Roxygen comments are shown below for those who read the scripts.
#
# @param param A vector of parameters in the FCVAR model.
# @param k The number of lags in the system.
# @param r The cointegrating rank.
# @param p The number of variables in the system.
# @param opt An S3 object of class \code{FCVAR_opt} that stores the chosen estimation options,
# generated from \code{FCVARoptions()}.
# @return \code{coeffs}, a list of coefficients of the FCVAR model.
# It has the same form as an element of the list \code{results}
# returned by \code{FCVARestn}.
#
# @examples
# opt <- FCVARoptions()
# opt$gridSearch   <- 0 # Disable grid search in optimization.
# opt$dbMin        <- c(0.01, 0.01) # Set lower bound for d,b.
# opt$dbMax        <- c(2.00, 2.00) # Set upper bound for d,b.
# opt$constrained  <- 0 # Impose restriction dbMax >= d >= b >= dbMin ? 1 <- yes, 0 <- no.
# x <- votingJNP2014[, c("lib", "ir_can", "un_can")]
# results <- FCVARestn(x, k = 2, r = 1, opt)
# params <- SEmat2vecU(coeffs = results$coeffs, k = 2, r = 1, p = 3, opt)
# coeffs <- SEvec2matU(param = params, k = 2, r = 1, p = 3, opt )
#
# params <- matrix(seq(25))
# coeffs <- SEvec2matU(param = params, k = 2, r = 1, p = 3, opt )
# params <- SEmat2vecU(coeffs = coeffs, k = 2, r = 1, p = 3, opt)
# @family FCVAR auxiliary functions
# @seealso \code{FCVARoptions} to set default estimation options.
# \code{FCVARestn} to estimate the FCVAR model and calculate
# standard errors of the estimates.
# \code{SEmat2vecU} is called by \code{FCVARhess} to convert the parameters
# from a vector into the coefficients after calculating the Hessian matrix.
# \code{SEmat2vecU} is a near inverse of \code{SEvec2matU},
# in the sense that \code{SEvec2matU} obtains only a
# subset of the parameters in \code{results$coeffs}.
#
SEvec2matU <- function(param, k, r, p, opt ) {

  # Create list for output.
  coeffs <- list(
    db = NULL,
    muHat = NULL,
    xiHat = NULL,
    alphaHat = NULL,
    GammaHat = NULL
  )

  if (opt$restrictDB) {
    coeffs$db <- matrix(c(param[1], param[1]), nrow = 1, ncol = 2)  # store d,b
    param <- param[2:length(param)]				# drop d,b from param
  } else {
    coeffs$db <- param[1:2]
    param <- param[3:length(param)]
  }


  if (opt$levelParam) {
    coeffs$muHat <- param[1:p]
    param <- param[(p+1):length(param)]
  }


  if (opt$unrConstant) {
    coeffs$xiHat <- matrix(param[1:p], nrow = p, ncol = 1)
    param <- param[(p+1):length(param)]
  }


  if (r > 0) {
    coeffs$alphaHat <- matrix( param[1:p*r], nrow = p, ncol = r)
    param <- param[(p*r+1):length(param)]
  } else {
    coeffs$alphaHat <- NULL
  }


  if (k > 0) {
    coeffs$GammaHat <- matrix( param[1 : length(param)], nrow = p, ncol = p*k)
  } else {
    coeffs$GammaHat <- NULL
  }


  return(coeffs)
}




# Calculate Restricted Estimates for the FCVAR Model
#
# \code{GetRestrictedParams} calculates restricted estimates of
# cointegration parameters and error variance in the FCVAR model.
# To calculate the optimum, it uses the switching algorithm
# of Boswijk and Doornik (2004, page 455) to optimize over free parameters
# \eqn{\psi} and \eqn{\phi} directly, combined with the line search proposed by
# Doornik (2016, working paper). We translate between  \eqn{(\psi, \phi)} and
# \eqn{(\alpha, \beta)} using the relation of \eqn{R_\alpha vec(\alpha) = 0} and
# \eqn{A\psi = vec(\alpha')}, and \eqn{R_\beta vec(\beta) = r_\beta} and
# \eqn{H\phi + h = vec(\beta)}. Note the transposes.
# Note that Roxygen comments are excluded to keep this function internal.
# However, the contents of Roxygen comments are shown below for those who read the scripts.
#
# @param beta0 The unrestricted estimate of \code{beta},
# a \eqn{p x r} matrix of cointegrating vectors
# returned from \code{FCVARestn} or \code{GetParams}.
# @param S00 A matrix of product moments,
# calculated from the output of \code{TransformData} in \code{GetParams}.
# @param S01 A matrix of product moments,
# calculated from the output of \code{TransformData} in \code{GetParams}.
# @param S11 A matrix of product moments,
# calculated from the output of \code{TransformData} in \code{GetParams}.
# @param cap_T The number of observations in the sample.
# @param p The number of variables in the system.
# @param opt An S3 object of class \code{FCVAR_opt} that stores the chosen estimation options,
# generated from \code{FCVARoptions()}.
#
# @return A list object \code{switched_mats} containing the restricted estimates,
# including the following parameters:
# \describe{
#   \item{\code{betaStar}}{A \eqn{p x r} matrix of cointegrating vectors.
#       The \eqn{r x 1} vector \eqn{\beta x_t} is the stationary cointegration relations.}
#   \item{\code{alphaHat}}{A \eqn{p x r} matrix of adjustment parameters.}
#   \item{\code{OmegaHat}}{A \eqn{p x p} covariance matrix of the error terms.}
# }
# @examples
# opt <- FCVARoptions()
# opt$gridSearch   <- 0 # Disable grid search in optimization.
# opt$dbMin        <- c(0.01, 0.01) # Set lower bound for d,b.
# opt$dbMax        <- c(2.00, 2.00) # Set upper bound for d,b.
# opt$constrained  <- 0 # Impose restriction dbMax >= d >= b >= dbMin ? 1 <- yes, 0 <- no.
# opt$R_Beta <- matrix(c(1, 0, 0), nrow = 1, ncol = 3)
# opt <- FCVARoptionUpdates(opt, p = 3, r = 1) # Need to update restriction matrices.
# betaStar <- matrix(c(-0.95335616, -0.07345676, -0.29277318), nrow = 3)
# S00 <- matrix(c(0.0302086527,  0.001308664,  0.0008200923,
#                 0.0013086640,  0.821417610, -0.1104617893,
#                 0.0008200923, -0.110461789,  0.0272861128), nrow = 3)
# S01 <- matrix(c(-0.0047314320, -0.04488533,  0.006336798,
#                 0.0026708007,  0.17463884, -0.069006455,
#                 -0.0003414163, -0.07110324,  0.022830494), nrow = 3, byrow = TRUE)
# S11 <- matrix(c( 0.061355941,  -0.4109969,  -0.007468716,
#                  -0.410996895,  70.6110313, -15.865097810,
#                  -0.007468716, -15.8650978,   3.992435799), nrow = 3)
# switched_mats <- GetRestrictedParams(betaStar, S00, S01, S11, cap_T = 316, p = 3, opt)
# @family FCVAR auxiliary functions
# @seealso \code{FCVARoptions} to set default estimation options.
# \code{FCVARestn} calls \code{GetParams} to estimate the FCVAR model,
# which in turn calls \code{GetRestrictedParams} if there are restrictions
# imposed on \code{alpha} or \code{beta}.
# @references Boswijk, H. P. and J. A. Doornik (2004).
# "Identifying, estimating and testing restricted cointegrated systems:
# An overview," Statistica Neerlandica 58, 440-465.
# @references Doornik, J. A. (2018).
# "Accelerated estimation of switching algorithms: the cointegrated
# VAR model and other applications,"
# Forthcoming in Scandinavian Journal of Statistics.
# @export
#
GetRestrictedParams <- function(beta0, S00, S01, S11, cap_T, p, opt) {


  r  <- ncol(beta0)
  p1 <- p + opt$rConstant

  # Restrictions on beta.
  if(is.null(opt$R_Beta)) {
    H <- diag(p1*r)
    h <- matrix(0, nrow = p1*r, ncol = 1)
  } else {
    H <- pracma::null(opt$R_Beta)
    h <- t(opt$R_Beta) %*%
      solve(opt$R_Beta %*% t(opt$R_Beta)) %*% opt$r_Beta
  }


  # Restrictions on alpha.
  #   We use the commutation matrix K_pr to transform vec(A) into vec(A'),
  #   see Magnus & Neudecker (1988, p. 47, eqn (1)).
  Ip  <- diag(p)
  Kpr <- matrix(kronecker(Ip, diag(r)), nrow = p*r, ncol = p*r)
  if(is.null(opt$R_Alpha)) {
    A <- Kpr %*% diag(p*r)
  } else {
    A <- pracma::null(opt$R_Alpha %*% solve(Kpr))
  }


  # Least squares estimator of Pi, used in calculations below.
  PiLS <- t((S01 %*% solve(S11)))
  vecPiLS <- matrix(PiLS, nrow = length(PiLS), ncol = 1)





  # Starting values for switching algorithm.
  betaStar <- beta0
  alphaHat <-  S01 %*% beta0 %*% solve(t(beta0) %*% S11 %*% beta0)
  OmegaHat <- S00 - S01 %*% betaStar %*% t(alphaHat) -
    alphaHat %*% t(betaStar) %*% t(S01) +
    alphaHat %*% t(betaStar) %*% S11 %*% betaStar %*% t(alphaHat)

  # Algorithm specifications.
  # iters   <- opt$UncFminOptions$MaxFunEvals
  # Tol     <- opt$UncFminOptions$TolFun
  iters   <- opt$unc_optim_control$maxit
  Tol     <- opt$unc_optim_control$reltol
  conv    <- 0
  i       <- 0

  # Tolerance for entering the line search epsilon_s in Doornik's paper.
  TolSearch <- 0.01

  # Line search parameters.
  lambda <- c(1, 1.2, 2, 4, 8)
  nS <- length(lambda)
  likeSearch <- matrix(NA, nrow = nS, ncol = 1)
  OmegaSearch <- array(0, dim = c(p, p,nS))


  # Get candidate values for entering the switching algorithm.
  vecPhi1 <- solve(t(H) %*%
                     kronecker(t(alphaHat) %*% solve(OmegaHat) %*% alphaHat, S11) %*%
                     H) %*%
    t(H) %*% (kronecker(t(alphaHat) %*% solve(OmegaHat), S11)) %*%
    (vecPiLS - kronecker(alphaHat, diag(p1)) %*% h)

  # Translate vecPhi to betaStar.
  vecB <- H %*% vecPhi1 + h
  betaStar <- matrix(vecB, nrow = p1, ncol = r)

  # Candidate value of vecPsi.
  vecPsi1 <- solve(t(A) %*%
                     kronecker(solve(OmegaHat), t(betaStar) %*% S11 %*% betaStar) %*%
                     A) %*%
    t(A) %*% (kronecker(solve(OmegaHat), t(betaStar) %*% S11)) %*% vecPiLS

  # Translate vecPsi to alphaHat.
  vecA <- A %*% vecPsi1 # This is vec(alpha')
  alphaHat <- matrix(solve(Kpr) %*% vecA, nrow = p, ncol = r)

  # Candidate values of piHat and OmegaHat.
  piHat1 <- alphaHat %*% t(betaStar)
  OmegaHat <- S00 - S01 %*% betaStar %*% t(alphaHat) -
    alphaHat %*% t(betaStar) %*% t(S01) +
    alphaHat %*% t(betaStar) %*% S11 %*% betaStar %*% t(alphaHat)

  # Calculate the likelihood.
  like1 <- - log(det(OmegaHat))

  while(i<=iters && !conv) {

    # Update values for convergence criteria.
    piHat0  <- piHat1
    like0   <- like1

    if(i == 1) {
      # Initialize candidate values.
      vecPhi0_c <- vecPhi1
      vecPsi0_c <- vecPsi1
    }


    #  ---- alpha update step ---- %
    # Update vecPsi.
    vecPsi1 <- solve(t(A) %*% kronecker(solve(OmegaHat), t(betaStar) %*%
                                     S11 %*% betaStar) %*% A) %*%
      t(A) %*% (kronecker(solve(OmegaHat), t(betaStar) %*% S11)) %*% vecPiLS

    # Translate vecPsi to alphaHat.
    vecA <- A %*% vecPsi1 # This is vec(alpha')
    alphaHat <- matrix(solve(Kpr) %*% vecA, nrow = p, ncol = r)

    #  ---- omega update step ---- %
    # Update OmegaHat.
    OmegaHat <- S00 - S01 %*% betaStar %*% t(alphaHat) -
      alphaHat %*% t(betaStar) %*% t(S01) +
      alphaHat %*% t(betaStar) %*% S11 %*% betaStar %*% t(alphaHat)

    #  ---- beta update step ---- %

    # Update vecPhi.
    vecPhi1 <- solve(t(H) %*%
                       kronecker(t(alphaHat) %*% solve(OmegaHat) %*% alphaHat, S11) %*%
                       H) %*%
      t(H) %*% (kronecker(t(alphaHat) %*% solve(OmegaHat), S11)) %*%
      (vecPiLS - kronecker(alphaHat, diag(p1)) %*% h)

    # Translate vecPhi to betaStar.
    vecB <- H %*% vecPhi1 + h
    betaStar <- matrix(vecB, nrow = p1, ncol = r)

    #  ---- pi and likelihood update  ---- %
    # Update estimate of piHat.
    piHat1 <- alphaHat %*% t(betaStar)

    # Update OmegaHat with new alpha and new beta.
    OmegaHat <- S00 - S01 %*% betaStar %*% t(alphaHat) -
      alphaHat %*% t(betaStar) %*% t(S01) +
      alphaHat %*% t(betaStar) %*% S11 %*% betaStar %*% t(alphaHat)

    # Calculate the likelihood.
    like1 <- - log(det(OmegaHat))

    if(i > 0) {

      # Calculate relative change in likelihood
      likeChange <- (like1 - like0) / (1 + abs(like0))

      # Check relative change and enter line search if below tolerance.
      if(likeChange < TolSearch && opt$LineSearch) {

        # Calculate changes in parameters.
        deltaPhi <- vecPhi1 - vecPhi0_c
        deltaPsi <- vecPsi1 - vecPsi0_c

        # Initialize parameter bins
        vecPhi2 <- matrix(NA, nrow = length(vecPhi1), ncol = nS)
        vecPsi2 <- matrix(NA, nrow = length(vecPsi1), ncol = nS)

        # Values already calculated for lambda = 1.
        vecPhi2[ ,1]       <- vecPhi1
        vecPsi2[ ,1]       <- vecPsi1
        likeSearch[1]      <- like1
        OmegaSearch[ , ,1] <- OmegaHat

        for (iL in 2:nS) {

          # New candidates for parameters based on line search.
          vecPhi2[ , iL] <- vecPhi0_c + lambda[iL]*deltaPhi
          vecPsi2[ , iL] <- vecPsi0_c + lambda[iL]*deltaPsi

          # Translate to alpha and beta
          vecA <- A %*% vecPsi2[ , iL]
          alphaHat <- matrix(solve(Kpr) %*% vecA, nrow = p, ncol = r)
          vecB <- H %*% vecPhi2[ , iL] + h
          betaStar <- matrix(vecB, nrow = p1, ncol = r)

          # Calculate and store OmegaHat
          OmegaSearch[ , , iL] <- S00 -
            S01 %*% betaStar %*% t(alphaHat) -

            alphaHat %*% t(betaStar) %*% t(S01) +
            alphaHat %*% t(betaStar) %*% S11 %*% betaStar %*% t(alphaHat)

          # Calculate and store log-likelihood
          likeSearch[iL] <- - log(det(OmegaSearch[ , , iL]))

        }

        # Update max likelihood and OmegaHat based on line search.
        iSearch <- which(likeSearch == max(likeSearch), arr.ind = TRUE)

        # If there are identical likelihoods, choose smallest
        # increment
        iSearch  <- min(iSearch)
        like1    <- likeSearch[iSearch]
        OmegaHat <- OmegaSearch[ , , iSearch]
        # Save old candidate parameter vectors for next iteration.
        vecPhi0_c  <- vecPhi1
        vecPsi0_c  <- vecPsi1

        # Update new candidate parameter vectors.
        vecPhi1  <- vecPhi2[ ,iSearch]
        vecPsi1  <- vecPsi2[ ,iSearch]
        # Update coefficients.
        vecA     <- A %*% vecPsi1
        alphaHat <- matrix(solve(Kpr) %*% vecA, nrow = p, ncol = r)
        vecB     <- H %*% vecPhi1 + h
        betaStar <- matrix(vecB, nrow = p1, ncol = r)
        # Update estimate of piHat.
        piHat1 <- alphaHat %*% t(betaStar)


      }


      # Calculate relative change in likelihood
      likeChange <- (like1 - like0) / (1 + abs(like0))

      # Calculate relative change in coefficients.
      piChange   <- max(abs(piHat1 - piHat0) / (1 + abs(piHat0)) )

      # Check convergence.
      if(abs(likeChange) <= Tol && piChange <= sqrt(Tol)) {
        conv <- 1
      }


    }

    i <- i+1
  }




  # Return a list of parameters from the switching algorithm.
  switched_mats <- list(
    betaStar = betaStar,
    alphaHat = alphaHat,
    OmegaHat = OmegaHat
  )

  return(switched_mats)
}
Any scripts or data that you put into this service are public.
FCVAR documentation built on May 5, 2022, 9:06 a.m.
rdrr.io home R language documentation Run R code online
CRAN packages Bioconductor packages R-Forge packages GitHub packages
Note that we can't provide technical support on individual packages. You should contact the package authors for that.
FCVAR
Estimation and Inference for the Fractionally Cointegrated VAR

R/FCVAR_aux.R
In FCVAR: Estimation and Inference for the Fractionally Cointegrated VAR

Defines functions GetRestrictedParams SEvec2matU SEmat2vecU FCVARhess GetFreeParams FracDiff Lbk GetResiduals TransformData FCVARlikeFull GetEstimates FCVARlike FCVARlikeMu find_local_max plot.FCVAR_grid FCVARlikeGrid GetParams FCVARsimBS FCVARsim

Documented in FCVARlikeGrid FCVARsim FCVARsimBS FracDiff plot.FCVAR_grid

Try the FCVAR package in your browser

R Package Documentation

Browse R Packages

We want your feedback!

FCVAR Estimation and Inference for the Fractionally Cointegrated VAR

R/FCVAR_aux.R In FCVAR: Estimation and Inference for the Fractionally Cointegrated VAR

Defines functions GetRestrictedParams SEvec2matU SEmat2vecU FCVARhess GetFreeParams FracDiff Lbk GetResiduals TransformData FCVARlikeFull GetEstimates FCVARlike FCVARlikeMu find_local_max plot.FCVAR_grid FCVARlikeGrid GetParams FCVARsimBS FCVARsim

Documented in FCVARlikeGrid FCVARsim FCVARsimBS FracDiff plot.FCVAR_grid

Try the FCVAR package in your browser

R Package Documentation

Browse R Packages

We want your feedback!

FCVAR
Estimation and Inference for the Fractionally Cointegrated VAR

R/FCVAR_aux.R
In FCVAR: Estimation and Inference for the Fractionally Cointegrated VAR
