Package {bayesianOU}

Type:	Package
Title:	Bayesian Nonlinear Ornstein-Uhlenbeck Models with Stochastic Volatility
Version:	0.2.0
Description:	Fits Bayesian nonlinear Ornstein-Uhlenbeck models with cubic drift, stochastic volatility, and Student-t innovations. The package implements hierarchical priors for sector-specific parameters and supports parallel MCMC sampling via 'Stan'. Model comparison is performed using Pareto Smoothed Importance Sampling Leave-One-Out (PSIS-LOO) cross-validation following Vehtari, Gelman, and Gabry (2017) <doi:10.1007/s11222-016-9696-4>. Prior specifications follow recommendations from Gelman (2006) <doi:10.1214/06-BA117A> for scale parameters.
License:	MIT + file LICENSE
Encoding:	UTF-8
Depends:	R (≥ 4.1.0)
Imports:	stats, graphics, utils, parallel
Suggests:	cmdstanr, rstan (≥ 2.21.0), loo (≥ 2.5.0), posterior, ggplot2, tidyr, openxlsx, knitr, rmarkdown, testthat (≥ 3.0.0)
Additional_repositories:	https://mc-stan.org/r-packages/
VignetteBuilder:	knitr
Config/testthat/edition:	3
URL:	https://github.com/IsadoreNabi/bayesianOU
BugReports:	https://github.com/IsadoreNabi/bayesianOU/issues
Config/roxygen2/version:	8.0.0
NeedsCompilation:	no
Packaged:	2026-06-18 18:57:36 UTC; josemgomezj
Author:	José Mauricio Gómez Julián [aut, cre]
Maintainer:	José Mauricio Gómez Julián <isadore.nabi@pm.me>
Repository:	CRAN
Date/Publication:	2026-06-18 19:30:02 UTC

bayesianOU: Bayesian Nonlinear Ornstein-Uhlenbeck Models

Description

Fits Bayesian nonlinear Ornstein-Uhlenbeck models with cubic drift, stochastic volatility (SV), and Student-t innovations. Sector-specific parameters use a non-centered hierarchical (partial-pooling) parameterization, and sampling runs via 'Stan' (with within-chain parallelism through reduce_sum).

Main Functions

• fit_ou_nonlinear_tmg: Fit the main OU model

• extract_posterior_summary: Extract posterior summaries

• validate_ou_fit: Validate model fit (MCMC + LOO + OOS)

• kappa_stability_evidence: Dynamic mean-reversion evidence

• compare_models_loo: Compare models via PSIS-LOO

Model Specification

The estimated one-step (Euler, dt = 1) mean equation is

\Delta Y_{t,s} = \kappa_s(\theta_s - Y_{t-1,s} + a_3 (Y_{t-1,s}-\theta_s)^3) + (\beta_{0,s} + \beta_1 zTMG_t) X_{t-1,s} + \gamma\, COM_{t-1,s} + \epsilon_{t,s}

with \epsilon \sim \mathrm{Student\text{-}t}(\nu, 0, \sigma_{t,s}) and \log\sigma^2_{t,s} a stationary AR(1). See the README for the full methodology (train/test design, priors, and the a_3<0 assumption).

Author(s)

Maintainer: José Mauricio Gómez Julián isadore.nabi@pm.me (ORCID)

Authors:

• José Mauricio Gómez Julián isadore.nabi@pm.me (ORCID)

See Also

Useful links:

• https://github.com/IsadoreNabi/bayesianOU

• Report bugs at https://github.com/IsadoreNabi/bayesianOU/issues

Align matrix columns to reference

Description

Reorders and filters columns of matrix B to match column names of A.

Usage

align_columns(A, B, verbose = FALSE)

Arguments

Value

Matrix B with columns reordered to match A

Build accounting block for TMG

Description

Creates accounting information for TMG transformations.

Usage

build_accounting_block(
  TMG_raw,
  zTMG_exo,
  zTMG_use,
  mu_tmg,
  sd_tmg,
  hard,
  sigma_delta
)

Arguments

Value

List with components:

tmg_byK

Back-transformed TMG used in model

tmg_exo

Back-transformed exogenous TMG

wedge_delta

Difference (zero if hard=TRUE)

sigma_delta_prior

Prior SD for wedge

note

Description of constraint type

Build beta(TMG_t) table by sector and time

Description

Constructs the time-varying beta matrix using posterior medians.

Usage

build_beta_tmg_table(fit, zTMG_use, summ = NULL)

Arguments

Details

This is a deterministic point reconstruction \beta_{s}(t) = \beta_{0,s} + \beta_1 \cdot zTMG_t evaluated at the posterior medians of beta0_s and beta1. It is NOT a sampled time-varying coefficient: the time variation comes only from zTMG_t.

Value

List with components:

beta_point

Matrix (T x S) of beta values

meta

List with description metadata

Check availability of Stan backend

Description

Verifies if cmdstanr or rstan is available for model fitting.

Usage

check_stan_backend(verbose = FALSE)

Arguments

Value

Character string: "cmdstanr", "rstan", or "none"

Compare models using PSIS-LOO

Description

Compares two fitted models using PSIS-LOO cross-validation.

Usage

compare_models_loo(results_new, results_base)

Arguments

Value

List with components:

loo_table

Comparison table from loo::loo_compare

deltaELPD

Numeric difference in ELPD

Examples

if (requireNamespace("loo", quietly = TRUE)) {
  # 1. Create mock 'loo' objects manually to avoid computation errors
  # Structure required by loo_compare: list with 'estimates' and class 'psis_loo'
  
  # Mock Model 1: ELPD = -100
  est1 <- matrix(c(-100, 5), nrow = 1, dimnames = list("elpd_loo", c("Estimate", "SE")))
  # Pointwise data (required for diff calculation). 10 observations.
  pw1 <- matrix(c(rep(-10, 10)), ncol = 1, dimnames = list(NULL, "elpd_loo"))
  
  loo_obj1 <- list(estimates = est1, pointwise = pw1)
  class(loo_obj1) <- c("psis_loo", "loo")
  
  # Mock Model 2: ELPD = -102 (worse)
  est2 <- matrix(c(-102, 5), nrow = 1, dimnames = list("elpd_loo", c("Estimate", "SE")))
  pw2 <- matrix(c(rep(-10.2, 10)), ncol = 1, dimnames = list(NULL, "elpd_loo"))
  
  loo_obj2 <- list(estimates = est2, pointwise = pw2)
  class(loo_obj2) <- c("psis_loo", "loo")
  
  # 2. Wrap in the structure expected by your package
  res_new <- list(diagnostics = list(loo = loo_obj1))
  res_base <- list(diagnostics = list(loo = loo_obj2))
  
  # 3. Compare (This will now run cleanly without warnings)
  cmp <- compare_models_loo(res_new, res_base)
  print(cmp)
}

Compute common factor from matrix

Description

Extracts first principal component and standardizes using training period.

Usage

compute_common_factor(Mz, T_train, use_train_loadings = TRUE, verbose = FALSE)

Arguments

Value

Numeric vector of factor scores (length T)

Count HMC divergences

Description

Extracts the number of divergent transitions from a fitted Stan model.

Usage

count_divergences(fit)

Arguments

Value

Integer. Number of divergent transitions (post-warmup).

Examples

# Create a "mock" CmdStanMCMC object for demonstration
# (This simulates a model with 0 divergences)
mock_fit <- structure(list(
  sampler_diagnostics = function() {
    # Return a 3D array: [iterations, chains, variables]
    # Variable 1 is usually accept_stat__, let's say var 2 is divergent__
    ar <- array(0, dim = c(100, 4, 6)) 
    dimnames(ar)[[3]] <- c("accept_stat__", "divergent__", "energy__", 
                           "n_leapfrog__", "stepsize__", "treedepth__")
    return(ar)
  }
), class = "CmdStanMCMC")

# Now the example can run without errors:
n_div <- count_divergences(mock_fit)
print(n_div)

Drift decomposition over grid

Description

Computes the cubic OU drift function over a grid of centered deviations z = Y - \theta_s. The evaluated function is \kappa_s(-z + a_{3,s} z^3), i.e. the drift expressed in the deviation coordinate, not in the original Y coordinate. To plot against Y, shift the grid by the corresponding theta_s.

Usage

drift_decomposition_grid(fit, summ, z_grid = seq(-2.5, 2.5, length.out = 101))

Arguments

Value

List with components:

z

The input z grid

drift

Matrix (length(z) x S) of drift values by sector

Evaluate out-of-sample forecast metrics

Description

Computes RMSE and MAE for multiple forecast horizons.

Usage

evaluate_oos(
  summ,
  Yz,
  Xz,
  zTMG,
  T_train,
  COM_ts,
  K_ts,
  com_in_mean = FALSE,
  horizons = c(1, 4, 8)
)

Arguments

Details

Two caveats matter for interpretation:

1. Conditional forecast. The recursion uses the realized future covariates X_{t-1} and zTMG_t at each step. It is therefore a forecast of Y conditional on the future paths of X and TMG, not an unconditional forecast. If those covariates are not known ex ante, treat these numbers as conditional/nowcasting metrics.

2. Plug-in medians. The path is propagated with posterior medians of the parameters and ignores parameter uncertainty, the stochastic volatility and the Student-t innovations. It is a point (mean-equation) forecast, not the full posterior predictive distribution. For genuinely out-of-sample numbers, fit with fit_window = "train".

Value

Named list with one element per horizon, each a list with h, RMSE, MAE and n_obs (number of pooled sector-by-origin errors; 0 / NA when the horizon exceeds the test window).

Examples

# 1. Generate dummy data for testing
T_obs <- 20
S <- 2
Yz <- matrix(rnorm(T_obs * S), nrow = T_obs, ncol = S)
Xz <- matrix(rnorm(T_obs * S), nrow = T_obs, ncol = S)
COM_ts <- matrix(abs(rnorm(T_obs * S)), nrow = T_obs, ncol = S)
K_ts <- matrix(abs(rnorm(T_obs * S)) + 1, nrow = T_obs, ncol = S)
zTMG <- rnorm(T_obs)

# 2. Create a dummy summary list (mimicking extract_posterior_summary)
summ <- list(
  theta_s = runif(S),
  kappa_s = runif(S),
  a3_s = runif(S),
  beta0_s = runif(S),
  beta1 = 0.5,
  gamma = 0.1
)

# 3. Run the function
metrics <- evaluate_oos(summ, Yz, Xz, zTMG, T_train = 15, 
                        COM_ts, K_ts, horizons = c(1, 2))
print(metrics)

Out-of-sample forecast metrics for the 2-level nested model

Description

Out-of-sample analogue of evaluate_oos for the nested model. The decisive difference from the single-level recursion is the attractor: in the single-level model the market price reverts to a constant level with a linear forcing by the exogenous production-price covariate X; in the nested model the market price \varphi reverts to the latent production price \Phi, which is itself propagated forward through its own Level-2 Ornstein-Uhlenbeck equation toward the G'-driven mean \mu_{s,t} = m_{0,s} + m_1 G'_t. Both states are advanced jointly; the market never sees the realized future X.

Usage

evaluate_oos_nested(
  meds,
  Yz,
  Phi_med,
  Gprime_z,
  zTMG,
  T_train,
  COM_ts,
  K_ts,
  kappa_cap,
  com_in_mean = FALSE,
  horizons = c(1, 4, 8),
  V_z = NULL
)

Arguments

Details

At each forecast origin t_0 the recursion is seeded with the realized standardized market price \varphi_{t_0-1} and the posterior-median latent production price \hat\Phi_{t_0-1}, then for h = 1,\dots steps:

\Phi_{t} = \Phi_{t-1} + \kappa^p_s(\mu_{s,t} - \Phi_{t-1}),\quad \mu_{s,t} = m_{0,s} + m_1 G'_t,

\varphi_{t} = \varphi_{t-1} + \kappa^m_s(-d_{t-1} + a_{3,s} d_{t-1}^3) + \gamma\, \mathrm{COM}_{t-1,s},\quad d_{t-1} = \varphi_{t-1} - \Phi_{t-1},

with \kappa^m_s = \texttt{kappa\_cap}\cdot \mathrm{logit}^{-1}(\tilde\kappa_s + \beta_1 zTMG_t) the bounded gravitation speed (the same logit link used in fitting). The market step uses \Phi_{t-1} (one-step lag), matching the Stan likelihood.

For the Level-3 value model (n_levels = 3), the Level-2 attractor also tracks the value index, \mu_{s,t} = m_{0,s} + m_1 G'_t + m_v V_{s,t} (production prices gravitate around values, Capital III ch. IX), recovered by supplying the standardized value anchor V_z and the coupling meds$m_v. Without them the recursion is the plain Level-2 mean.

Two caveats, identical in spirit to evaluate_oos:

1. Conditional forecast. The recursion still consumes the realized future G'_t and zTMG_t (the macro drivers), so it is a forecast conditional on those aggregate paths, not on the sectoral production prices (which are now latent and forecast endogenously).

2. Plug-in medians. States are propagated with posterior medians and ignore parameter uncertainty, the latent-\Phi innovation, the stochastic volatility and the Student-t tails. It is a point (mean-equation) forecast, not the full posterior predictive distribution. Fit with fit_window = "train" for genuinely held-out numbers.

Value

Named list with one element per horizon, each a list with h, RMSE, MAE and n_obs (pooled sector-by-origin errors; 0 / NA when the horizon exceeds the test window). Errors are in standardized \varphi units, directly comparable with evaluate_oos.

See Also

evaluate_oos (single-level), fit_ou_nested.

Export posterior draws to CSV

Description

Saves posterior draws for selected parameters to a CSV file.

Usage

export_draws_csv(fit, params, path, verbose = FALSE)

Arguments

Value

Path to the created file, invisibly.

Export model comparison to Excel

Description

Creates an Excel workbook with model comparison results, parameter summaries, and fit information.

Usage

export_model_comparison(
  results_new,
  results_base,
  path = file.path(tempdir(), "model_comparison.xlsx"),
  verbose = FALSE
)

Arguments

Details

Creates three worksheets:

• loo_comparison: PSIS-LOO comparison table

• param_summary: Sector-specific parameter estimates

• fit_info: Model configuration and diagnostics

Value

TRUE invisibly on success.

Examples

if (requireNamespace("openxlsx", quietly = TRUE) &&
    requireNamespace("loo", quietly = TRUE)) {
  
  # 1. Create mock results objects
  # Mock Model New
  res_new <- list(
    factor_ou = list(
      kappa_s = c(0.5, 0.6), a3_s = c(-0.1, -0.2), beta0_s = c(1, 2),
      gamma = 0.05, model = "TestModel", beta1 = 0.3, nu = 4,
      factor_ou_info = list(T_train = 50, com_in_mean = TRUE)
    ),
    diagnostics = list(
      divergences = 0,
      # Mock LOO object
      loo = list(
         estimates = matrix(c(-100, 2), 1, 2, dimnames=list("elpd_loo", c("Estimate","SE"))),
         pointwise = matrix(rep(-2, 50), ncol=1)
      )
    )
  )
  class(res_new$diagnostics$loo) <- c("psis_loo", "loo")
  
  # Mock Model Base
  res_base <- list(
    diagnostics = list(
      loo = list(
         estimates = matrix(c(-110, 2), 1, 2, dimnames=list("elpd_loo", c("Estimate","SE"))),
         pointwise = matrix(rep(-2.2, 50), ncol=1)
      )
    )
  )
  class(res_base$diagnostics$loo) <- c("psis_loo", "loo")

  # 2. Define a safe temporary path
  out_path <- file.path(tempdir(), "comparison.xlsx")
  
  # 3. Run export (This writes to tempdir, allowed by CRAN)
  try({
    export_model_comparison(res_new, res_base, path = out_path)
    # unlink(out_path) # Cleanup
  })
}

Export model summary to text file

Description

Creates a plain text summary of model results.

Usage

export_summary_txt(fit_res, path, verbose = FALSE)

Arguments

Value

Path to the created file, invisibly.

Mean-reversion (dynamic stability) evidence for kappa parameters

Description

Computes 95 percent credible intervals for each kappa_s (mean reversion speed) and the posterior probability that every sector reverts.

Usage

extract_convergence_evidence(fit_res, verbose = TRUE)

Arguments

Value

List with components:

kappa_ic95

Matrix (S x 3) with columns q2.5, median, q97.5

stable

Logical indicating if all kappa CIs fall in (0,1)

convergence

Deprecated alias of stable (kept for back-compat)

prob_stable

Posterior probability of joint mean reversion

prob_convergence

Deprecated alias of prob_stable

Important - this is NOT MCMC convergence

Despite the historical name, this function does NOT assess sampler convergence. It evaluates a dynamic property of the estimated process: whether the mean-reversion speed lies in a range consistent with a stable, reverting (gravitating) system. For sampler convergence use R-hat, ESS and divergences (see validate_ou_fit). The threshold kappa < 1 is a (conservative) monotone-reversion criterion; under the Euler discretization the linear map is stable for 0 < kappa < 2. The alias kappa_stability_evidence is preferred.

Examples

# 1. Create a mock fit object with kappa draws
# kappa_tilde is log(kappa), so we use log(0.5) roughly -0.69
n_draws <- 100
S <- 2
kappa_tilde_draws <- matrix(rnorm(n_draws * S, mean = -0.7, sd = 0.1), 
                            nrow = n_draws, ncol = S)
colnames(kappa_tilde_draws) <- c("kappa_tilde[1]", "kappa_tilde[2]")

mock_fit <- structure(list(
  draws = function(vars, format="matrix") {
    if (vars == "kappa_tilde") return(kappa_tilde_draws)
    return(NULL)
  }
), class = "CmdStanMCMC")

# 2. Wrap in the results list structure
results_mock <- list(
  factor_ou = list(
    stan_fit = mock_fit
  )
)

# 3. Extract evidence
conv <- extract_convergence_evidence(results_mock)
print(conv$prob_convergence)

Latent Level-2 mean trajectory mu_t = m0_s + m1 * G'_t

Description

Reconstructs the time-varying Level-2 mean toward which the latent production price \Phi reverts, \mu_{s,t} = m_{0,s} + m_1 G'_t, with credible bands obtained by evaluating the expression over the posterior draws of mu_const (m_{0,s}) and m1. This is the structural hook to the (not-yet-formalized) Level 3 of values: \mu is the slow attractor of the production price, here driven by the aggregate profit rate G'.

Usage

extract_mu_trajectory(fit, Gprime_z, sector_names = NULL, V_z = NULL)

Arguments

Value

List with median, q2.5, q97.5 (each a T x S matrix) and Gprime_z (the driver), or NULL if the Level-2 parameters are not present in the fit.

See Also

fit_ou_nested.

Extract posterior summary from fitted model

Description

Extracts median point estimates and credible intervals for all model parameters from a fitted Stan model.

Usage

extract_posterior_summary(fit)

Arguments

Value

List with components:

beta1

Median of global TMG effect

beta0_s

Vector of sector-specific intercepts

kappa_s

Vector of mean reversion speeds

a3_s

Vector of cubic drift coefficients

theta_s

Vector of equilibrium levels

rho_s

Vector of SV persistence parameters

alpha_s

Vector of SV level parameters

sigma_eta_s

Vector of SV volatility parameters

nu

Degrees of freedom for Student-t

gamma

COM effect in mean

rhat

R-hat convergence diagnostics

ess

Effective sample sizes

Examples

# 1. Create a mock CmdStanMCMC object
# We simulate a posterior distribution for 2 sectors
S <- 2
n_draws <- 100

# Helper to generate random draws
mock_draws <- function(name, n_cols=1) {
  m <- matrix(rnorm(n_draws * n_cols), nrow = n_draws, ncol = n_cols)
  if (n_cols > 1) {
    colnames(m) <- paste0(name, "[", 1:n_cols, "]")
  } else {
    colnames(m) <- name
  }
  as.data.frame(m)
}

# Combine draws into one data frame
df_draws <- cbind(
  mock_draws("beta1", 1),
  mock_draws("beta0_s", S),
  mock_draws("kappa_tilde", S), # Note: function expects log-scale kappa
  mock_draws("a3_tilde", S),    # Note: function expects log-scale a3
  mock_draws("theta_s", S),
  mock_draws("rho_s", S),
  mock_draws("alpha_s", S),
  mock_draws("sigma_eta_s", S),
  mock_draws("nu_tilde", 1),
  mock_draws("gamma", 1)
)

# Mock fit object
mock_fit <- structure(list(
  draws = function(vars, format="df") {
     # Simple regex matching for the mock
     if (length(vars) == 1) {
       # Check if it's a scalar or vector parameter request
       if (vars %in% names(df_draws)) return(df_draws[vars])
       # Pattern match for vectors like "beta0_s" -> "beta0_s[1]", "beta0_s[2]"
       cols <- grep(paste0("^", vars, "\\["), names(df_draws), value = TRUE)
       if (length(cols) > 0) return(df_draws[cols])
     }
     return(df_draws) 
  },
  summary = function() {
    data.frame(
      variable = names(df_draws),
      rhat = rep(1.0, ncol(df_draws)),
      ess_bulk = rep(400, ncol(df_draws))
    )
  }
), class = "CmdStanMCMC")

# 2. Run extraction
summ <- extract_posterior_summary(mock_fit)
print(summ$kappa_s)

Fit the unified nonlinear OU model: single-level or 2-level nested

Description

Single engine behind the package. n_levels = 1 reproduces the legacy single-level model exactly (market price reverts to a constant level with a linear forcing by the production price and a TMG interaction); n_levels = 2 fits the nested model in which the market price \varphi reverts to a latent production price \Phi (Level 1), and \Phi has its own Ornstein-Uhlenbeck dynamics reverting to a G'-driven mean (Level 2). The constructed production-price index enters as a noisy measurement (anchor) of the latent \Phi.

Usage

fit_ou_nested(
  Y,
  X,
  TMG,
  COM,
  CAPITAL_TOTAL,
  n_levels = 1L,
  Gprime = NULL,
  k_cost = NULL,
  V_value = NULL,
  level_spec = NULL,
  theta_separation = c("soft", "hard"),
  k_uncertainty = c("meas", "recon"),
  sigma_phi_meas_fixed = NULL,
  kappa_cap = 2,
  results_robust = list(),
  priors = list(),
  com_in_mean = TRUE,
  train_frac = 0.7,
  fit_window = c("train", "full"),
  chains = 6,
  iter = 12000,
  warmup = 6000,
  thin = 1,
  cores = max(1, parallel::detectCores() - 1),
  threads_per_chain = 2,
  hard_sum_zero = TRUE,
  orthogonalize_tmg = TRUE,
  factor_from = c("X", "Y"),
  use_train_loadings = TRUE,
  adapt_delta = 0.97,
  max_treedepth = 12,
  seed = 1234,
  init = NULL,
  moment_match = NULL,
  verbose = FALSE
)

Arguments

Details

fit_ou_nonlinear_tmg is a thin backward-compatible wrapper around this engine with n_levels = 1; this function is the single source of truth.

Value

A list. For n_levels = 1 the structure matches the legacy single-level fit, with the inner factor_ou carrying the historical S3 class "ou_nonlinear_tmg". For n_levels = 2 the top-level object carries the S3 class "ou_nested_2level" (so print, summary and plot dispatch on it) and adds factor_ou$level2 (production reversion speed kappa_p, market speed at zTMG=0 kappa_m_base, mean intercept mu_const, slope m1, innovation scale sigma_p, measurement scale sigma_phi_meas, and the Level-2 richness quantities that are present only when their switch is on: the cubic restoring coefficient a3_p (l2_cubic), the Student-t degrees of freedom nu_p (l2_studentt) and the time-varying SV scale path sigma_p_t = \exp(h_p/2) (l2_sv); see decision D-IMPL-10); phi_latent (the latent \Phi trajectory with credible bands); mu_path (the G'-driven Level-2 mean trajectory \mu_{s,t} = m_{0,s} + m_1 G'_t with bands); a separation block with the posterior evidence for \kappa^m > \kappa^p; and diagnostics$oos, the 2-level out-of-sample metrics in which the market reverts to the latent \Phi propagated forward (see evaluate_oos_nested).

See Also

fit_ou_nonlinear_tmg for the single-level wrapper and fit_ou_nested_mi for the multiple-imputation driver that couples disaggregation draws to the 2-level model via Rubin's rule.

Examples

T_obs <- 24; S <- 2
Y <- matrix(rnorm(T_obs * S), T_obs, S)
X <- matrix(rnorm(T_obs * S), T_obs, S)
TMG <- rnorm(T_obs)
COM <- matrix(runif(T_obs * S), T_obs, S)
K <- matrix(runif(T_obs * S, 100, 1000), T_obs, S)
Gp <- rnorm(T_obs)
if (requireNamespace("cmdstanr", quietly = TRUE)) {
  try(fit_ou_nested(Y, X, TMG, COM, K, n_levels = 2, Gprime = Gp,
                    chains = 1, iter = 100, warmup = 50), silent = TRUE)
}

Fit the nested OU model under multiple imputation of the market price

Description

Couples the disaggregation posterior to the OU model by multiple imputation (Rubin's rule). The disaggregation step produces M posterior draws of the sector-level market price \varphi (each a T x S matrix); this driver refits the OU model once per imputation and combines the analyses, propagating the disaggregation uncertainty into the OU posterior. Two complementary summaries are returned: Rubin's moment-pooled estimates (within + between imputation variance) and the mixture posterior over the imputations (the "complete-draws" pooling used for joint probabilities such as the time-scale separation).

Usage

fit_ou_nested_mi(
  phi_draws,
  X,
  TMG,
  COM,
  CAPITAL_TOTAL,
  Gprime = NULL,
  V_value = NULL,
  n_levels = 2L,
  M = 25,
  pool_params = NULL,
  keep_fits = FALSE,
  keep_kappa_mixture = FALSE,
  kappa_mixture_draws = 4000L,
  seed = 1234,
  verbose = FALSE,
  checkpoint_dir = NULL,
  ...
)

Arguments

Details

Unified external multiple imputation (D-IMPL-13.5). Beyond the market price, X (the production-price anchor \Phi), V_value (the value anchor), CAPITAL_TOTAL (K) and COM may also vary by imputation. Each of these four accepts either a fixed matrix (the legacy behaviour: used verbatim for every imputation, bit-identical results) or a per-imputation set – a list of D matrices or a 3-D array [T, S, D] – paired one-to-one with phi_draws. This lets a single external imputation generator propagate ALL the construction uncertainty (disaggregation of the market, the imputed split of v and p that feeds K and V, and the reconstructed \Phi) through the convergent K-DET / Variant-A engine, pooling by Rubin's rule. The aggregates TMG and Gprime stay fixed (they are observed, preserved by the aggregate-preserving renormalization of the generator). The engine is untouched (R-side only): n \in \{1, 2\} and the all-fixed path remain bit-identical.

Value

An object of class "ou_nested_mi": a list with

rubin

Data frame of Rubin-pooled estimates per parameter (estimate, total SD, 95% CI, Barnard-Rubin df, fraction of missing information).

phi_latent_pooled

(n_levels >= 2) Rubin-pooled latent \Phi path: mean, sd, q2.5, q97.5 as T x S matrices.

separation_pooled

(n_levels >= 2) mixture-posterior evidence for \kappa^m > \kappa^p: prob_sep_by_sector (per-sector P(\kappa^m > \kappa^p), the sectoral distribution to report – at n_levels = 3 the joint probability typically collapses because \kappa^p is re-identified once the value term absorbs the production mean, see D-IMPL-10.1, which is NOT "no separation"), prob_sep_joint, and the per-sector reversion summaries kappa_m_median / kappa_p_median and half-lives (\ln 2)/\kappa (halflife_*_median, halflife_*_q2.5, halflife_*_q97.5) from the mixture posterior. Since \kappa is positive by construction, P(\kappa>0)=1 is uninformative; the half-life (and whether it exceeds the observed span – the low-kappa trap) is the substantive per-sector quantity.

kappa_mixture

(only when keep_kappa_mixture = TRUE) a thinned mixture posterior of kappa_m_base and kappa_p (draws x S), stacked across imputations, for sectoral half-life distributions and low-kappa-trap probabilities at a chosen horizon.

per_imputation

List of compact per-imputation summaries or full fits if keep_fits = TRUE. Each summary carries the chain-aware convergence over the pooled parameters (rhat_max_pooled, ess_bulk_min_pooled, ess_tail_min_pooled) – the binding diagnostic for Rubin's pooling – plus the global rhat_max, rhat_share, ess_bulk_min (over all parameters, latent states included; reference only) and divergences.

config

The imputation/pooling configuration.

References

Rubin DB (1987) Multiple Imputation for Nonresponse in Surveys. Barnard J, Rubin DB (1999) Small-sample degrees of freedom with multiple imputation. Biometrika 86(4):948-955.

See Also

fit_ou_nested.

Fit Bayesian nonlinear OU model with TMG effect and SV (single-level)

Description

Backward-compatible single-level entry point. This is a thin wrapper around fit_ou_nested with n_levels = 1: the nested engine is the single source of truth, and n_levels = 1 reproduces the legacy single-level model exactly (market price reverts to a constant level with a linear forcing by the production price X and a TMG interaction). The returned object keeps the historical structure and S3 class "ou_nonlinear_tmg" so existing extractors, plots and validators work unchanged.

Usage

fit_ou_nonlinear_tmg(
  results_robust,
  Y,
  X,
  TMG,
  COM,
  CAPITAL_TOTAL,
  model = c("base"),
  priors = list(),
  com_in_mean = TRUE,
  train_frac = 0.7,
  fit_window = c("train", "full"),
  chains = 6,
  iter = 12000,
  warmup = 6000,
  thin = 1,
  cores = max(1, parallel::detectCores() - 1),
  threads_per_chain = 2,
  hard_sum_zero = TRUE,
  orthogonalize_tmg = TRUE,
  factor_from = c("X", "Y"),
  use_train_loadings = TRUE,
  adapt_delta = 0.97,
  max_treedepth = 12,
  seed = 1234,
  init = NULL,
  moment_match = NULL,
  verbose = FALSE
)

Arguments

Details

The model uses a non-centered hierarchical (partial-pooling) parameterization for the sector-specific parameters theta_s, kappa_s, a3_s and beta0_s: each is built as hyper_intercept + hyper_slope * COM_s + hyper_sd * z, with z ~ N(0,1). The training window is floor(T * train_frac). All data standardization, the COM weighting, and (by default) the common factor loadings are computed from the training window only. The likelihood window is controlled by fit_window to keep train/test coherent.

The cubic drift enforces a3 < 0 (a restoring force that strengthens with the deviation); this is a stability assumption, not an estimated result, and it precludes detecting locally expansive (self-amplifying) regimes. The Euler discretization uses dt = 1, so the discrete persistence of the linear part is 1 - kappa; interpret half-lives accordingly.

Value

List containing factor_ou (model results and parameter estimates), beta_tmg, sv, nonlinear, accounting and diagnostics (MCMC diagnostics, PSIS-LOO and OOS metrics). See fit_ou_nested for the underlying engine and the 2-level model.

Examples

# 1. Prepare dummy data
T_obs <- 20
S_sectors <- 2
Y <- matrix(rnorm(T_obs * S_sectors), nrow = T_obs, ncol = S_sectors)
X <- matrix(rnorm(T_obs * S_sectors), nrow = T_obs, ncol = S_sectors)
TMG <- rnorm(T_obs)
COM <- matrix(runif(T_obs * S_sectors), nrow = T_obs, ncol = S_sectors)
K <- matrix(runif(T_obs * S_sectors, 100, 1000), nrow = T_obs, ncol = S_sectors)

# 2. Run model (conditional on Stan backend availability)
if (requireNamespace("cmdstanr", quietly = TRUE) ||
    requireNamespace("rstan", quietly = TRUE)) {
  try({
    results <- fit_ou_nonlinear_tmg(
      results_robust = list(),
      Y = Y, X = X, TMG = TMG, COM = COM, CAPITAL_TOTAL = K,
      chains = 1, iter = 100, warmup = 50,
      verbose = FALSE
    )
  }, silent = TRUE)
}

Mean-reversion (dynamic stability) evidence for kappa parameters

Description

Preferred alias of extract_convergence_evidence. See that function for details. The name avoids conflating dynamic mean reversion with MCMC sampler convergence.

Usage

kappa_stability_evidence(fit_res, verbose = TRUE)

Arguments

Value

See extract_convergence_evidence.

Bridge a methods-enabled cmdstan model to the geometry engine

Description

Decoupled core (the improvement over the gdpar bridge, which required a fitted object of a specific class and recompiled from its source): takes a CmdStanModel ALREADY compiled with compile_model_methods = TRUE together with its Stan data, samples one iteration to expose the standalone log_prob/grad_log_prob/hessian methods, reads the unconstrained dimension and a posterior-mean warm-start, and returns an ou_geom_target plus a reference position. Nothing in the regular fit path is touched.

Usage

ou_geom_bridge(
  model,
  stan_data,
  dim = NULL,
  reference = NULL,
  hessian = TRUE,
  methods_seed = 1L
)

Arguments

Value

A list with target (an ou_geom_target), dim, reference, fit (the one-iteration cmdstan fit with methods) and has_hessian.

Static Hamiltonian Monte Carlo with a pluggable geometry

Description

Fixed step-size / fixed trajectory-length HMC with a Metropolis correction over a pluggable metric. The default Euclidean metric is textbook HMC; a Riemannian metric reuses the same loop with the generalised implicit leapfrog of Girolami & Calderhead (2011). The metric is a preconditioner only, so the returned draws are exact regardless of which metric is chosen.

Usage

ou_geom_hmc(
  target,
  metric = NULL,
  epsilon = 0.1,
  L = 20L,
  n_iter = 1000L,
  n_warmup = 500L,
  init = NULL,
  seed = NULL,
  fp_tol = 1e-09,
  fp_max = 100L
)

Arguments

Value

A list of class ou_geom_hmc with draws, accept_rate, n_divergent, energy, ebfmi, epsilon, L and metric_type.

Euclidean (constant) metric for the geometry engine

Description

A position-independent mass matrix: the identity (diagonal) by default, or a supplied symmetric positive-definite matrix / positive variance vector (dense preconditioner; the remedy for a straight anisotropic canyon).

Usage

ou_geom_metric_euclidean(dim = NULL, M = NULL)

Arguments

Value

A ou_geom_metric (position-independent).

Riemannian (position-dependent) metric for the geometry engine

Description

A position-dependent mass M(\theta) adapting the sampler's local distance to the curvature of the log-posterior – the remedy for a funnel (variable curvature). curvature = "softabs" maps the eigenvalues of the Hessian of -\log\pi through \lambda\coth(\alpha\lambda) (Betancourt 2013), flooring nearly flat directions; it needs only the Hessian (from the target or finite-differenced). curvature = "fisher" uses a supplied expected-Fisher function (positive-definite by construction).

Usage

ou_geom_metric_riemannian(
  target,
  curvature = c("softabs", "fisher"),
  fisher = NULL,
  dfisher = NULL,
  alpha = 1e+06,
  floor = 1e-08,
  fd_step = 1e-04
)

Arguments

Value

A ou_geom_metric (position-dependent).

Build a geometry-engine sampling target

Description

Wrap an unconstrained log-density (and its gradient, optionally its Hessian) into the target object ou_geom_hmc integrates. Accepts either a cmdstan fit/model compiled with compile_model_methods = TRUE (whose $log_prob / $grad_log_prob / $hessian act on the unconstrained scale) or explicit R closures.

Usage

ou_geom_target(
  object = NULL,
  log_prob = NULL,
  grad_log_prob = NULL,
  dim = NULL,
  hessian = NULL,
  param_names = NULL,
  data = NULL
)

Arguments

Value

An object of class ou_geom_target.

Named per-level richness configurations for the comparative experiment

Description

Convenience builder for the level_spec argument of fit_ou_nested. Each level has four richness switches: cubic (cubic drift), sv (stochastic volatility), student_t (Student-t innovations) and hierarchy (cross-sector partial pooling). The switches act only in the 2-level mode.

Usage

ou_level_spec(
  config = c("canonical", "n1_full_n2_lean", "both_full", "both_lean", "n1_lean",
    "n1_lean_n2_lean")
)

Arguments

Details

The comparative experiment of the design contrasts these configurations by parameter recovery and out-of-sample / leave-cluster-out predictive performance, reporting which one the data support (no configuration is claimed uniquely correct):

"canonical" (alias "n1_full_n2_lean")

Level 1 full (cubic + SV + Student-t + hierarchy), Level 2 lean (linear Gaussian OU with a hierarchical mean). The default.

"both_full"

Both levels full: the latent production price also gets cubic restoring, stochastic volatility and Student-t innovations.

"both_lean"

Both levels lean: linear Gaussian OU at each level, hierarchy retained.

"n1_lean" (alias "n1_lean_n2_lean")

Level 1 lean, Level 2 lean.

The fifth experiment arm, the single-level model, is obtained directly with fit_ou_nested(..., n_levels = 1) (it takes no level_spec).

Value

A level_spec list with level1 and level2 entries, each a list of the four logical switches, ready to pass to fit_ou_nested.

See Also

fit_ou_nested.

Examples

ou_level_spec("both_full")
ou_level_spec("both_lean")

Stan code for the unified nonlinear OU model

Description

Returns the complete Stan code for the unified Ornstein-Uhlenbeck model: a single file (inst/stan/ou_nested.stan) covering the single-level mode (n_levels = 1, cubic drift, stochastic volatility, Student-t innovations, non-centered hierarchical priors and a TMG interaction) and the 2-level nested mode (n_levels = 2, market price reverting to a latent production price with its own G'-driven OU). The code is read from the canonical file (single source of truth); this function does not embed a duplicate copy.

Usage

ou_nested_stan_code()

Value

Character string containing the Stan model code.

Examples

code <- ou_nested_stan_code()
cat(substr(code, 1, 300))

Plot a 2-level nested OU fit

Description

Plot a 2-level nested OU fit

Usage

## S3 method for class 'ou_nested_2level'
plot(x, type = c("phi", "mu", "sv_p", "separation"), sector = 1, ...)

Arguments

Value

NULL invisibly; draws as a side effect.

Plot a multiple-imputation nested OU fit

Description

Plot a multiple-imputation nested OU fit

Usage

## S3 method for class 'ou_nested_mi'
plot(x, type = c("phi", "fmi"), sector = 1, ...)

Arguments

Value

NULL invisibly; draws as a side effect.

Plot beta(TMG_t) evolution by sector

Description

Creates a line plot showing the evolution of time-varying beta coefficients for selected sectors.

Usage

plot_beta_tmg(results, sectors = NULL)

Arguments

Value

A ggplot2 object if ggplot2 is available, otherwise NULL with a base R plot produced as side effect.

Examples

# 1. Create mock data (T x S matrix)
T_obs <- 50
S <- 3
beta_mat <- matrix(rnorm(T_obs * S), nrow = T_obs, ncol = S)
colnames(beta_mat) <- paste0("Sector_", 1:S)

# 2. Wrap in list structure expected by function
results_mock <- list(
  beta_tmg = list(
    beta_point = beta_mat
  )
)

# 3. Plot
plot_beta_tmg(results_mock)

Plot cubic OU drift curves

Description

Displays the drift function for selected sectors showing mean reversion with cubic nonlinearity.

Usage

plot_drift_curves(results, sectors = NULL)

Arguments

Value

NULL invisibly. Produces a base R plot as side effect.

Examples

# 1. Create mock data
# z: vector of state deviations
z_seq <- seq(-3, 3, length.out = 100)
# drift: matrix (rows=z, cols=sectors)
drift_mat <- cbind(
  -0.5 * z_seq - 0.1 * z_seq^3, # Sector 1
  -0.8 * z_seq - 0.05 * z_seq^3 # Sector 2
)

# 2. Wrap in list structure
results_mock <- list(
  nonlinear = list(
    drift_decomp = list(
      z = z_seq,
      drift = drift_mat
    )
  )
)

# 3. Plot
plot_drift_curves(results_mock)

Plot posterior distributions of key parameters

Description

Creates density plots for selected parameters.

Usage

plot_posterior_densities(fit, params = c("beta1", "nu"), verbose = FALSE)

Arguments

Value

NULL invisibly. Produces plots as side effect.

Plot stochastic volatility evolution

Description

Displays the estimated volatility path for a selected sector.

Usage

plot_sv_evolution(results, sector = 1)

Arguments

Value

NULL invisibly. Produces a base R plot as side effect.

Examples

# 1. Create mock data (Volatility must be positive)
T_obs <- 50
sigma_mat <- matrix(exp(rnorm(T_obs * 2)), nrow = T_obs, ncol = 2)

# 2. Wrap in list structure
results_mock <- list(
  sv = list(
    h_summary = list(
      sigma_t = sigma_mat
    )
  )
)

# 3. Plot sector 1
plot_sv_evolution(results_mock, sector = 1)

Print a 2-level nested OU fit

Description

Compact one-screen overview: data dimensions, MCMC health, the time-scale separation evidence \kappa^m > \kappa^p, and the median ranges of the market and production reversion speeds.

Usage

## S3 method for class 'ou_nested_2level'
print(x, ...)

Arguments

Value

x, invisibly.

Print a multiple-imputation nested OU fit

Description

Print a multiple-imputation nested OU fit

Usage

## S3 method for class 'ou_nested_mi'
print(x, ...)

Arguments

Value

x, invisibly.

Summarize stochastic volatility sigmas

Description

Extracts median volatility paths from the SV component.

Usage

summarize_sv_sigmas(fit)

Arguments

Value

List with component:

sigma_t

Matrix (T x S) of median volatilities

Summarize a 2-level nested OU fit

Description

Assembles the Level-2 parameter table (production reversion speed, market speed at zTMG=0, mean intercept, innovation scale, measurement scale), the separation evidence, the MCMC diagnostics and the out-of-sample metrics.

Usage

## S3 method for class 'ou_nested_2level'
summary(object, ...)

## S3 method for class 'summary.ou_nested_2level'
print(x, ...)

Arguments

Value

An object of class summary.ou_nested_2level.

Summarize a multiple-imputation nested OU fit

Description

Summarize a multiple-imputation nested OU fit

Usage

## S3 method for class 'ou_nested_mi'
summary(object, ...)

## S3 method for class 'summary.ou_nested_mi'
print(x, ...)

Arguments

Value

An object of class summary.ou_nested_mi.

Validate OU model fit

Description

Prints diagnostic summaries and hypothesis tests for a fitted model.

Usage

validate_ou_fit(fit_res, verbose = TRUE)

Arguments

Value

Invisibly returns the diagnostics list

Examples

# Create a dummy results list that mimics the output of fit_ou_nonlinear_tmg
dummy_results <- list(
  diagnostics = list(
    rhat = c(alpha = 1.01, beta = 1.00),
    ess = c(alpha = 400, beta = 350),
    loo = list(estimates = matrix(c(1, 0.1), ncol=2, 
               dimnames=list("elpd_loo", c("Estimate", "SE")))),
    oos = list(h1 = list(RMSE = 0.5))
  ),
  factor_ou = list(beta1 = 0.3),
  nonlinear = list(a3 = -0.5),
  sv = list(rho_s = 0.2)
)

# Run validation on the dummy object
validate_ou_fit(dummy_results)

Verbose message helper

Description

Prints a message only if verbose is TRUE.

Usage

vmsg(msg, verbose)

Arguments

Value

NULL invisibly

Z-score standardization using training period statistics

Description

Standardizes a matrix using mean and standard deviation computed from the training period only.

Usage

zscore_train(M, T_train, eps = 1e-08)

Arguments

Value

A list with components:

Mz

Standardized matrix of same dimensions as M

mu

Vector of training means for each column

sd

Vector of training standard deviations for each column

Examples

M <- matrix(rnorm(100), nrow = 20, ncol = 5)
result <- zscore_train(M, T_train = 14)
str(result)