R/03_simulate_paths.R

In weatherMRJD: Weather Analysis and Markov Regime Switching Jump Diffusion Models

#' Simulate Temperature Paths via Non-Homogeneous Jump-Diffusion
#'
#' Generates discretized Monte Carlo sample trajectories of daily average temperatures.
#'
#' @param n_paths Integer count of simulation paths.
#' @param days Horizon length in days (default 365).
#' @param kappa Mean-reversion speed.
#' @param sigma_0 Baseline volatility.
#' @param lambda_0 Baseline jump intensity.
#' @param lambda_1 Amplitude of dynamic seasonal jump intensity.
#' @param mu_jump Mean jump magnitude.
#' @param sigma_jump Volatility of jump magnitude.
#' @param baseline_temp Baseline temperature parameter A.
#' @param amplitude_temp Seasonal temperature amplitude.
#' @return An (n_paths x days) matrix of simulated daily temperatures.
#' @examples
#' paths <- simulate_weather_paths(n_paths = 50, days = 30)
#' dim(paths)
#' @export
simulate_weather_paths <- function(n_paths = 5000, days = 365, kappa = 0.08, 
                                   sigma_0 = 2.5, lambda_0 = 0.05, lambda_1 = 0.03,
                                   mu_jump = -3.0, sigma_jump = 1.5,
                                   baseline_temp = 14.5, amplitude_temp = 11.2) {
  t <- seq(0, days - 1, length.out = days)
  S_t <- baseline_temp + amplitude_temp * sin(2.0 * pi * (t - 80.0) / 365.0)
  sigma_t  <- sigma_0 * (1.0 + 0.3 * cos(2.0 * pi * t / 365.0))
  lambda_t <- lambda_0 + lambda_1 * cos(2.0 * pi * (t - 180.0) / 365.0)
  
  X <- matrix(0, nrow = n_paths, ncol = days)
  T_paths <- matrix(0, nrow = n_paths, ncol = days)
  T_paths[, 1] <- S_t[1]
  
  for (day in 2:days) {
    Z <- stats::rnorm(n_paths)
    p_jump <- 1.0 - exp(-lambda_t[day] * 1.0)
    jumps_occur <- stats::runif(n_paths) < p_jump
    jump_sizes  <- stats::rnorm(n_paths, mean = mu_jump, sd = sigma_jump) * jumps_occur
    
    dX <- -kappa * X[, day - 1] + sigma_t[day] * Z + jump_sizes
    X[, day] <- X[, day - 1] + dX
    T_paths[, day] <- S_t[day] + X[, day]
  }
  
  return(T_paths)
}