PREDICT INDONESIA 2025 INFLATION

Even worse, the precise result of a computation to which floating point truncation error applies is not necessarily invariant across systems. Operating system differences, compiler settings, and a host of other factors can influence the outcome. The details don’t matter for this post, and frankly I don’t understand all of them myself. For all I know the colour of the laptop case might be relevant, or the name of the programmer’s cat. Weirdness abounds once your calculations start to run up against the limits of floating point precision.

Just for the sake of argument then, let’s imagine that during the course of some fancy simulation, you and I compute a covariance matrix on different machines. It’s supposed to be the same covariance matrix, but thanks to the weirdness of floating point your machine computes cov1 and mine computes cov2. The differences are very small, but they’re large enough that – SOMEHOW, IN DEFIANCE OF ALL THE LAWS OF GOD AND MAN – this happens:

dsf

transformer <- function(sigma, method) {
  
  # extract necessary quantities from eigendecomposition 
  if (method %in% c("eigen-mass", "eigen-mvtnorm")) {
    eig <- eigen(sigma)
    Q <- eig$vectors             # matrix of eigenvectors
    L <- diag(sqrt(eig$values))  # diagonal matrix of sqrt eigenvalues
  }
  
  # compute the transformation matrix A for eigendecomposition
  if (method == "eigen-mass") A <- L %*% t(Q)           # MASS-style
  if (method == "eigen-mvtnorm") A <- Q %*% L %*% t(Q)  # mvtnorm-style
  
  # compute the transformation matrix A for cholesky with pivoting
  if (method == "chol-mvtnorm") {
    U <- chol(sigma, pivot = TRUE) # upper triangular matrix 
    ord <- order(attr(U, "pivot")) # reordering required due to pivot
    A <- U[, ord]
  }
  
  return(A)
}

These are of course not the only ways to do it but this post is already giving me body horror nightmares so I refuse to implement the matrix square root method. Not going to happen.

For our purposes it is useful to split it into two parts. First, we write a transformer() function that returns a transformation matrix that we can use when sampling from multivariate normal with covariance matrix :