Simulation DGPs

# Replication script for the optional e, shifts, coef_noise and coef_a axes of
# sim_psy1() and for sim_blan(type = "rotermann_wilfling") in R/sim.R.
devtools::load_all("exuber", quiet = TRUE)

cat("=== sim_psy1(e = ...): custom innovations replace rnorm(n-1, sd=sigma) ===\n")
set.seed(42)
y_plain <- sim_psy1(50, seed = 42)
y_zero_e <- sim_psy1(50, seed = 42, e = rep(0, 49))
cat("y_zero_e finite:", all(is.finite(y_zero_e)), "\n")
cat("y_zero_e != y_plain (e overrides rnorm):", !isTRUE(all.equal(y_zero_e, y_plain)), "\n\n")

cat("=== sim_psy1(shifts = ...): one-period deterministic bump ===\n")
y_shift <- sim_psy1(50, seed = 42, shifts = list(date = 10, size = 100))
cat("y_shift[10] - y_plain[10] (expect 100):", y_shift[10] - y_plain[10], "\n")
cat("y_shift[1:9] == y_plain[1:9] (no effect before shift date):",
    isTRUE(all.equal(y_shift[1:9], y_plain[1:9])), "\n\n")

cat("=== sim_blan(type = 'rotermann_wilfling') ===\n")
set.seed(7)
n_b <- 5; pi_b <- 0.7; delta_b <- 0.984; rw_sigma <- 0.05; b0 <- 0.1
theta_b <- rbinom(n_b - 1, 1, pi_b)
u_b <- rlnorm(n_b - 1, meanlog = -rw_sigma ^ 2 / 2, sdlog = rw_sigma)
cat("theta:", theta_b, "\n")
cat("u:", paste(sprintf("%.8f", u_b), collapse = ","), "\n")
b <- b0
for (i in 1:(n_b - 1)) {
  b[i + 1] <- if (theta_b[i] == 1) {
    b[i] * u_b[i] / delta_b
  } else {
    ((1 - pi_b * delta_b) / (1 - pi_b)) * b[i] * u_b[i]
  }
}
cat("b:", paste(sprintf("%.8f", b), collapse = ","), "\n\n")

cat("=== sim_blan(type = 'rotermann_wilfling') via the real function ===\n")
set.seed(7)
b_real <- sim_blan(n_b, pi = pi_b, type = "rotermann_wilfling", delta = delta_b,
                    rw_sigma = rw_sigma, b0 = b0)
cat("b (real fn call):", paste(sprintf("%.8f", b_real), collapse = ","), "\n")
cat("matches hand-traced b:", isTRUE(all.equal(as.numeric(b_real), b)), "\n")

# All numbers above are reproduced (the deterministic/formula parts
# bit-for-bit, the RNG-driven parts structurally) in
# docs/replication/simulation-dgps/sim_psy1_axes_and_blan_rw_validation.py
# and in pyexuber/tests/test_sim.py.

"""Python counterpart of sim_psy1_axes_and_blan_rw_validation.R --
cross-checks pyexuber's sim_psy1() e/shifts/coef_noise/coef_a optional
axes and sim_blan(type="rotermann_wilfling") (exuber.sim) against R's
R/sim.R equivalents.

Run standalone: uv run --project pyexuber python
docs/replication/simulation-dgps/sim_psy1_axes_and_blan_rw_validation.py

sim_psy1's e=/shifts= structural behavior (overrides innovations, adds
an isolated one-period bump) is checked directly against the real
function. sim_blan(type="rotermann_wilfling")'s recursion formula is
checked bit-for-bit via a fake RNG feeding R's exact theta/u draws to
the real sim_blan() function.
"""

import numpy as np

from exuber.sim import sim_blan, sim_psy1


def check_sim_psy1_e() -> None:
    y_plain = sim_psy1(50, seed=42)
    y_zero_e = sim_psy1(50, seed=42, e=np.zeros(49))
    assert np.all(np.isfinite(y_zero_e))
    assert not np.allclose(y_zero_e, y_plain)  # e overrides the default rnorm draws
    print("sim_psy1(e=...): overrides innovations, matches R's structural behavior.")


def check_sim_psy1_shifts() -> None:
    y_plain = sim_psy1(50, seed=42)
    y_shift = sim_psy1(50, seed=42, shifts={"date": [10], "size": [100.0]})
    assert y_shift[9] - y_plain[9] == 100.0  # R: y_shift[10] - y_plain[10] == 100
    assert np.array_equal(y_shift[:9], y_plain[:9])
    print("sim_psy1(shifts=...): exact +100 one-period bump, no effect before it, matches R.")


def check_sim_blan_rotermann_wilfling() -> None:
    # set.seed(7); theta <- rbinom(4,1,0.7); u <- rlnorm(4, meanlog=-.05^2/2, sdlog=.05) in R.
    theta = np.array([0, 1, 1, 1])
    u = np.array([0.96467442, 0.97837265, 0.95143523, 0.95254875])

    class _FakeRWRNG:
        def binomial(self, n_trials, p, size=None):
            return theta

        def lognormal(self, mean=0.0, sigma=1.0, size=None):
            return u

    orig = np.random.default_rng
    np.random.default_rng = lambda seed=None: _FakeRWRNG()
    try:
        b = sim_blan(
            5, pi=0.7, type="rotermann_wilfling", delta=0.984, rw_sigma=0.05, b0=0.1, seed=7
        )
    finally:
        np.random.default_rng = orig

    expected = [0.10000000, 0.10006889, 0.09949661, 0.09620385, 0.09312891]
    np.testing.assert_allclose(b, expected, atol=1e-6)
    print("sim_blan(type='rotermann_wilfling'): recursion matches R bit-for-bit.")


if __name__ == "__main__":
    check_sim_psy1_e()
    check_sim_psy1_shifts()
    check_sim_blan_rotermann_wilfling()
    print("All checks passed.")

# Replication script for the innovation generators in R/sim.R: sim_vol_break(),
# sim_vol_garch() (TGARCH through gamma > 0), sim_vol_cir(), sim_vol_sv(),
# sim_fi() and sim_innov().
#
# All of them are deterministic transformations of Gaussian noise, so each is
# checked by fixing the underlying rnorm() draws (set.seed(1); rnorm(10)) and
# tracing the recursion by hand. pyexuber/tests/test_sim.py uses the same
# technique for the regime-switching branches of sim_psy1().
devtools::load_all("exuber", quiet = TRUE)

set.seed(1)
eps10 <- rnorm(10)
cat("eps10:", paste(sprintf("%.8f", eps10), collapse = ","), "\n\n")

cat("=== sim_vol_break(n=10, tau=0.5, ratio=3, sigma=2) ===\n")
n <- 10; tau <- 0.5; ratio <- 3; sigma <- 2
sd_t <- sigma * ifelse(seq_len(n) > floor(tau * n), ratio, 1)
z_break <- sd_t * eps10
cat("sd_t:", sd_t, "\n")
cat("z:", paste(sprintf("%.8f", z_break), collapse = ","), "\n\n")

cat("=== sim_vol_garch(n=5, omega=0.1, alpha=0.1, beta=0.8, gamma=0) ===\n")
omega <- 0.1; alpha <- 0.1; beta <- 0.8
eps5 <- eps10[1:5]
h <- numeric(5); z <- numeric(5); h_prev <- 0; z_prev <- 0
for (t in 1:5) {
  h[t] <- omega + alpha * z_prev ^ 2 + beta * h_prev
  z[t] <- sqrt(h[t]) * eps5[t]
  h_prev <- h[t]; z_prev <- z[t]
}
cat("h:", paste(sprintf("%.8f", h), collapse = ","), "\n")
cat("z:", paste(sprintf("%.8f", z), collapse = ","), "\n\n")

cat("=== sim_vol_garch(..., gamma=0.2) [TGARCH] ===\n")
gamma2 <- 0.2
h2 <- numeric(5); z2 <- numeric(5); h_prev <- 0; z_prev <- 0
for (t in 1:5) {
  h2[t] <- omega + alpha * z_prev ^ 2 + beta * h_prev + gamma2 * z_prev ^ 2 * (z_prev < 0)
  z2[t] <- sqrt(h2[t]) * eps5[t]
  h_prev <- h2[t]; z_prev <- z2[t]
}
cat("h:", paste(sprintf("%.8f", h2), collapse = ","), "\n")
cat("z:", paste(sprintf("%.8f", z2), collapse = ","), "\n\n")

cat("=== sim_vol_cir(n=5, kappa=0.03, theta=0.25, xi=0.1) ===\n")
kappa <- 0.03; theta <- 0.25; xi <- 0.1; sigma0_sq <- theta
dt <- 1 / 5
db <- eps10[1:4] * sqrt(dt)
mult <- eps10[6:10]
sig2 <- numeric(5); sig2[1] <- sigma0_sq
for (i in 2:5) {
  prev <- max(sig2[i - 1], 0)
  sig2[i] <- max(prev + kappa * (theta - prev) * dt + xi * sqrt(prev) * db[i - 1], 0)
}
cat("sig2:", paste(sprintf("%.8f", sig2), collapse = ","), "\n")
cat("z:", paste(sprintf("%.8f", sqrt(sig2) * mult), collapse = ","), "\n\n")

cat("=== sim_vol_sv(n=5, phi=0.98, tau=0.1) ===\n")
phi <- 0.98; tau_sv <- 0.1
eta <- eps10[1:4] * tau_sv
mult_sv <- eps10[6:10]
log_sig2 <- numeric(5); log_sig2[1] <- 0
for (i in 2:5) log_sig2[i] <- phi * log_sig2[i - 1] + eta[i - 1]
cat("log_sig2:", paste(sprintf("%.8f", log_sig2), collapse = ","), "\n")
cat("z:", paste(sprintf("%.8f", exp(log_sig2 / 2) * mult_sv), collapse = ","), "\n\n")

cat("=== sim_fi(): psi recursion (d=0.2, m=5) + hand convolution (n=3) ===\n")
d <- 0.2; m <- 5
psi <- numeric(m + 1); psi[1] <- 1
for (j in 2:(m + 1)) psi[j] <- psi[j - 1] * (j - 2 + d) / (j - 1)
cat("psi:", paste(sprintf("%.8f", psi), collapse = ","), "\n")
epsfi <- eps10[1:8]
u <- vapply(1:3, function(k) sum(psi * rev(epsfi[k:(k + m)])), numeric(1))
cat("u:", paste(sprintf("%.8f", u), collapse = ","), "\n\n")

cat("=== sim_innov(dist='normal', sigma=2) ===\n")
cat("z:", paste(sprintf("%.8f", eps10 * 2), collapse = ","), "\n\n")

cat("=== sim_innov: t/skew_t closed-form constants (df=5) ===\n")
df <- 5
cat("t rescale 1/sqrt(df/(df-2)):", sprintf("%.10f", 1 / sqrt(df / (df - 2))), "\n")
xi <- -0.75
delta <- xi / sqrt(1 + xi ^ 2)
e_abs_t0 <- (2 * sqrt(df) / ((df - 1) * beta(df / 2, 0.5))) / sqrt(df / (df - 2))
mean_raw <- delta * e_abs_t0
sd_raw <- sqrt(max(1 - delta ^ 2 * e_abs_t0 ^ 2, .Machine$double.eps))
cat("delta:", sprintf("%.10f", delta), " e_abs_t0:", sprintf("%.10f", e_abs_t0),
    " mean_raw:", sprintf("%.10f", mean_raw), " sd_raw:", sprintf("%.10f", sd_raw), "\n")

# All numbers above are reproduced bit-for-bit in
# docs/replication/simulation-dgps/sim_vol_innovations_validation.py and
# in pyexuber/tests/test_sim.py (same reference values, duplicated as
# literals there per docs/replication/README.md's convention).

"""Python counterpart of sim_vol_innovations_validation.R -- cross-checks
pyexuber's sim_vol_break()/sim_vol_garch()/sim_vol_cir()/sim_vol_sv()/
sim_fi()/sim_innov() (exuber.sim) against R's R/sim.R equivalents.

Run standalone: uv run --project pyexuber python
docs/replication/simulation-dgps/sim_vol_innovations_validation.py

These are otherwise-deterministic transformations of Gaussian noise, so
each is checked by feeding the *actual* exuber.sim functions a fixed
pool of draws (R: set.seed(1); rnorm(10)) via a fake np.random.Generator
replacement (_FakeRNG) that returns pool slices in the same call order R
draws them in -- the algorithm gets checked bit-for-bit, not the RNG
bit-stream (see exuber.sim's module docstring on that distinction).
"""

import math

import numpy as np

from exuber.sim import (
    _beta_fn,
    sim_fi,
    sim_innov,
    sim_vol_break,
    sim_vol_cir,
    sim_vol_garch,
    sim_vol_sv,
)

# set.seed(1); rnorm(10) in R.
EPS10 = np.array(
    [-0.62645381, 0.18364332, -0.83562861, 1.59528080, 0.32950777,
     -0.82046838, 0.48742905, 0.73832471, 0.57578135, -0.30538839]
)


class _FakeRNG:
    """Feeds a fixed pool to standard_normal()/normal() in call order --
    see exuber's own tests/test_sim.py, which uses the identical class."""

    def __init__(self, pool):
        self.pool = np.asarray(pool, dtype=float)
        self.pos = 0

    def _take(self, size):
        if size is None:
            v = self.pool[self.pos]
            self.pos += 1
            return v
        out = self.pool[self.pos : self.pos + size]
        self.pos += size
        return out

    def standard_normal(self, size=None):
        return self._take(size)

    def normal(self, loc=0.0, scale=1.0, size=None):
        if size is None and isinstance(scale, np.ndarray):
            size = scale.shape[0]
        return loc + scale * self._take(size)


def _with_fake_rng(pool, fn, *args, **kwargs):
    orig = np.random.default_rng
    np.random.default_rng = lambda seed=None: _FakeRNG(pool)
    try:
        return fn(*args, **kwargs)
    finally:
        np.random.default_rng = orig


def check_sim_vol_break() -> None:
    y = _with_fake_rng(EPS10, sim_vol_break, 10, tau=0.5, ratio=3, sigma=2, seed=1)
    expected = [-1.25290762, 0.36728665, -1.67125722, 3.19056160, 0.65901554,
                -4.92281030, 2.92457431, 4.42994823, 3.45468811, -1.83233032]
    np.testing.assert_allclose(y, expected, atol=1e-6)
    print("sim_vol_break(): matches R.")


def check_sim_vol_garch() -> None:
    y = _with_fake_rng(EPS10[:5], sim_vol_garch, 5, omega=0.1, alpha=0.1, beta=0.8, gamma=0.0, seed=1)
    expected = [-0.19810209, 0.07875803, -0.41593815, 0.89607126, 0.21675026]
    np.testing.assert_allclose(y, expected, atol=1e-6)

    y2 = _with_fake_rng(EPS10[:5], sim_vol_garch, 5, omega=0.1, alpha=0.1, beta=0.8, gamma=0.2, seed=1)
    expected2 = [-0.19810209, 0.08042096, -0.42119779, 0.95247029, 0.22731252]
    np.testing.assert_allclose(y2, expected2, atol=1e-6)
    print("sim_vol_garch(): GARCH and TGARCH (gamma>0) match R.")


def check_sim_vol_cir() -> None:
    pool = np.concatenate([EPS10[:4], EPS10[5:10]])
    y = _with_fake_rng(pool, sim_vol_cir, 5, kappa=0.03, theta=0.25, xi=0.1, seed=1)
    expected = [-0.41023419, 0.23678823, 0.36175334, 0.27117724, -0.15439035]
    np.testing.assert_allclose(y, expected, atol=1e-6)
    print("sim_vol_cir(): matches R.")


def check_sim_vol_sv() -> None:
    pool = np.concatenate([EPS10[:4], EPS10[5:10]])
    y = _with_fake_rng(pool, sim_vol_sv, 5, phi=0.98, tau=0.1, seed=1)
    expected = [-0.82046838, 0.47239810, 0.72260999, 0.54069901, -0.31098370]
    np.testing.assert_allclose(y, expected, atol=1e-6)
    print("sim_vol_sv(): matches R.")


def check_sim_innov_normal() -> None:
    y = _with_fake_rng(EPS10, sim_innov, 10, dist="normal", sigma=2, seed=1)
    expected = [-1.25290762, 0.36728665, -1.67125722, 3.19056160, 0.65901554,
                -1.64093677, 0.97485810, 1.47664941, 1.15156270, -0.61077677]
    np.testing.assert_allclose(y, expected, atol=1e-6)
    print("sim_innov(dist='normal'): matches R.")


def check_sim_innov_t_skew_t_constants() -> None:
    df = 5
    assert 1 / math.sqrt(df / (df - 2)) == 0.7745966692 or abs(
        1 / math.sqrt(df / (df - 2)) - 0.7745966692
    ) < 1e-9

    xi = -0.75
    delta = xi / math.sqrt(1 + xi**2)
    e_abs_t0 = (2 * math.sqrt(df) / ((df - 1) * _beta_fn(df / 2, 0.5))) / math.sqrt(df / (df - 2))
    mean_raw = delta * e_abs_t0
    sd_raw = math.sqrt(max(1 - delta**2 * e_abs_t0**2, np.finfo(float).eps))
    assert abs(delta - (-0.6)) < 1e-9
    assert abs(e_abs_t0 - 0.7351051939) < 1e-9
    assert abs(mean_raw - (-0.4410631163)) < 1e-9
    assert abs(sd_raw - 0.8974760874) < 1e-9
    print("sim_innov(dist='t'/'skew_t'): closed-form rescaling constants match R.")


def check_sim_fi_psi_and_convolution() -> None:
    d, m = 0.2, 5
    psi = np.empty(m + 1)
    psi[0] = 1.0
    for j in range(1, m + 1):
        psi[j] = psi[j - 1] * (j - 1 + d) / j
    np.testing.assert_allclose(psi, [1.0, 0.2, 0.12, 0.088, 0.0704, 0.059136], atol=1e-6)

    eps = EPS10[:8]  # n=3, m=5 -> n+m=8 innovations, matching sim_fi's own convolution
    u = np.empty(3)
    for k in range(3):
        window = eps[k : k + m + 1]
        u[k] = np.dot(psi, window[::-1])
    np.testing.assert_allclose(u, [-0.66078593, 0.45529270, 0.82924303], atol=1e-6)
    print("sim_fi(): psi recursion and convolution formula match R.")

    # smoke test the real function (m = max(500, 5n) there, too big to fake-RNG by hand)
    y = sim_fi(200, d=0.2, sigma=1.0, seed=1)
    assert y.shape == (200,) and np.all(np.isfinite(y))
    print("sim_fi(): real call (n=200) runs and returns finite values.")


if __name__ == "__main__":
    check_sim_vol_break()
    check_sim_vol_garch()
    check_sim_vol_cir()
    check_sim_vol_sv()
    check_sim_innov_normal()
    check_sim_innov_t_skew_t_constants()
    check_sim_fi_psi_and_convolution()
    print("All checks passed.")