diff --git a/core/docs/plans/M3-sims-hub.md b/core/docs/plans/M3-sims-hub.md index aaf475a..2430b60 100644 --- a/core/docs/plans/M3-sims-hub.md +++ b/core/docs/plans/M3-sims-hub.md @@ -15,10 +15,11 @@ DESIGN-FIRST · ABSENT. Role C3; implementation C1. Mathematical foundations C4 established); specific model parameters C1. ## 3. Language & location -TBD · `src/economy/sims/`. Each sim type uses the runtime suited to its math: **Fortran** -(M3d, M3e — dense numerical PDE/SDE), **Prolog** (M3b, M3f — game-theoretic equilibria), -**R** (M3a — statistical inference), **Solidity** (M3c — on-chain precision), **Zig** -(M3g — tick-level latency). A query facade accessible to Traders. +**Tcl** · `src/economy/sims/`. The hub is a syntax-agnostic coordinator: Tcl manages +lifecycle, tick-advancement, and query routing for sub-sims in their native runtimes via +stdin/stdout JSON — **Fortran** (M3d, M3e), **Prolog** (M3b, M3f), **R** (M3a), +**Solidity** (M3c), **Zig** (M3g). Tcl imposes no type system or paradigm on the +sub-processes it orchestrates. ## 4. Does / does-not - **Does:** tick-advance continuously at **90:1** (1 wall-second = 90 simulated seconds) diff --git a/core/src/economy/sims/hub.tcl b/core/src/economy/sims/hub.tcl new file mode 100644 index 0000000..893a3f7 --- /dev/null +++ b/core/src/economy/sims/hub.tcl @@ -0,0 +1,284 @@ +#!/usr/bin/env tclsh +# +# M3 Sims Hub — lifecycle manager and query facade for M3a–M3g sims. +# Syntax-agnostic coordinator: speaks to R, Prolog, Solidity, Fortran, Zig +# sub-processes via stdin/stdout JSON. + +package require Tcl 8.6 + +namespace eval ::sims { + + variable SIM_TYPES { + statistical {lang R dir statistical} + sociological {lang Prolog dir sociological} + amm_liquidity {lang Solidity dir amm} + mev_adversarial {lang Fortran dir mev} + tokenomics_macro {lang Fortran dir tokenomics} + consensus_staking {lang Prolog dir consensus} + market_microstructure {lang Zig dir microstructure} + } + + variable SIM_PROCS + array set SIM_PROCS {} + + variable BASE_SPEED 90 + variable TICK_INTERVAL_MS 1000 + + # BoundedPrediction constructor — L2: every output has explicit bounds + proc bounded_prediction {value lower upper confidence horizon sim_type} { + if {$lower > $value || $value > $upper} { + error "invariant violation: lower <= value <= upper required\ + (got $lower <= $value <= $upper)" + } + if {$confidence < 0.0 || $confidence > 10.0} { + error "confidence must be in \[0.00, 10.00\], got $confidence" + } + dict create \ + value $value \ + lower_bound $lower \ + upper_bound $upper \ + confidence $confidence \ + time_horizon $horizon \ + sim_type $sim_type \ + timestamp [clock milliseconds] + } + + proc format_confidence {conf} { + format "%.2f/10.00" $conf + } + + proc format_gain {lower value upper horizon} { + format "%.2f - %.2f - %.2f / 10.00 gain over next %s" \ + $lower $value $upper $horizon + } + + # Prediction query structure + proc prediction_query {pred_type horizon {params {}}} { + dict create \ + prediction_type $pred_type \ + time_horizon $horizon \ + params $params \ + timestamp [clock milliseconds] + } + + # Sim status — L1: sims are always running + proc sim_status {sim_type} { + variable SIM_PROCS + if {[info exists SIM_PROCS($sim_type)]} { + set info $SIM_PROCS($sim_type) + dict create \ + running true \ + sim_type $sim_type \ + last_calibration [dict get $info last_cal] \ + data_freshness [expr {[clock milliseconds] - [dict get $info last_cal]}] \ + tick_count [dict get $info ticks] + } else { + dict create \ + running false \ + sim_type $sim_type + } + } + + # Query a running sim — L4: read-only, never mutates state + proc query {sim_type query_dict} { + variable SIM_PROCS + variable SIM_TYPES + + if {![dict exists $SIM_TYPES $sim_type]} { + error "unknown sim type: $sim_type\ + (valid: [dict keys $SIM_TYPES])" + } + + if {![info exists SIM_PROCS($sim_type)]} { + error "sim $sim_type is not running" + } + + set proc_info $SIM_PROCS($sim_type) + set chan [dict get $proc_info channel] + + set query_json [dict_to_json $query_dict] + puts $chan $query_json + flush $chan + + set response [gets $chan] + set result [json_to_dict $response] + + set pred [dict get $result prediction] + set bp [bounded_prediction \ + [dict get $pred value] \ + [dict get $pred lower_bound] \ + [dict get $pred upper_bound] \ + [dict get $pred confidence] \ + [dict get $pred time_horizon] \ + $sim_type] + + return $bp + } + + # Calibrate a sim with fresh data from M2 — L4: only M2 writes + proc calibrate {sim_type feed_data} { + variable SIM_PROCS + if {![info exists SIM_PROCS($sim_type)]} { + error "sim $sim_type is not running — cannot calibrate" + } + + set proc_info $SIM_PROCS($sim_type) + set chan [dict get $proc_info channel] + + set cal_msg [dict create \ + type calibrate \ + data $feed_data] + puts $chan [dict_to_json $cal_msg] + flush $chan + + dict set SIM_PROCS($sim_type) last_cal [clock milliseconds] + } + + # Launch a sim subprocess — L5: each sim is independent + proc launch_sim {sim_type} { + variable SIM_TYPES + variable SIM_PROCS + + if {![dict exists $SIM_TYPES $sim_type]} { + error "unknown sim type: $sim_type" + } + + set spec [dict get $SIM_TYPES $sim_type] + set lang [dict get $spec lang] + set dir [dict get $spec dir] + + set cmd [resolve_launcher $lang $dir] + + set chan [open "| $cmd" r+] + fconfigure $chan -buffering line -blocking 1 + + set SIM_PROCS($sim_type) [dict create \ + channel $chan \ + lang $lang \ + dir $dir \ + pid [pid $chan] \ + ticks 0 \ + last_cal [clock milliseconds] \ + started [clock milliseconds]] + + return [sim_status $sim_type] + } + + # Stop a sim — L5: failure in one doesn't cascade + proc stop_sim {sim_type} { + variable SIM_PROCS + if {[info exists SIM_PROCS($sim_type)]} { + set chan [dict get $SIM_PROCS($sim_type) channel] + catch {puts $chan {{"type":"shutdown"}}} + catch {close $chan} + unset SIM_PROCS($sim_type) + } + } + + # Resolve the launch command for a sim's language + proc resolve_launcher {lang dir} { + set base [file dirname [info script]] + switch -- $lang { + R { return "Rscript --vanilla ${base}/${dir}/main.R" } + Prolog { return "swipl -q -f ${base}/${dir}/main.pl" } + Solidity { return "node ${base}/${dir}/runner.js" } + Fortran { return "${base}/${dir}/sim" } + Zig { return "${base}/${dir}/sim" } + default { error "no launcher for language: $lang" } + } + } + + # Tick all running sims — L3: all horizons concurrent, 90:1 + proc tick_all {} { + variable SIM_PROCS + variable BASE_SPEED + + set tick_msg [dict create \ + type tick \ + sim_seconds $BASE_SPEED \ + wall_ms 1000] + + set tick_json [dict_to_json $tick_msg] + + foreach sim_type [array names SIM_PROCS] { + set chan [dict get $SIM_PROCS($sim_type) channel] + if {[catch { + puts $chan $tick_json + flush $chan + dict incr SIM_PROCS($sim_type) ticks + } err]} { + puts stderr "sim $sim_type tick failed: $err" + } + } + } + + # Minimal JSON serialization for Tcl dicts + proc dict_to_json {d} { + set pairs {} + dict for {k v} $d { + if {[string is double -strict $v]} { + lappend pairs "\"$k\":$v" + } elseif {[string is boolean -strict $v]} { + lappend pairs "\"$k\":[expr {$v ? "true" : "false"}]" + } elseif {[string index $v 0] eq "\{" || [string index $v 0] eq "\["} { + lappend pairs "\"$k\":$v" + } else { + lappend pairs "\"$k\":\"[string map {\" \\\" \\ \\\\} $v]\"" + } + } + return "\{[join $pairs ,]\}" + } + + proc json_to_dict {json} { + set json [string trim $json "\{\}"] + set d [dict create] + foreach pair [split $json ,] { + if {[regexp {"([^"]+)"\s*:\s*(.*)} $pair -> k v]} { + set v [string trim $v] + set v [string trim $v "\""] + dict set d $k $v + } + } + return $d + } + + # Main loop — L1: sims are always running, L3: tick-advanced + proc run_loop {} { + variable TICK_INTERVAL_MS + while {1} { + tick_all + after $TICK_INTERVAL_MS + } + } +} + +# Self-test when run directly +if {[info script] eq $::argv0} { + puts "M3 Sims Hub — Tcl [info patchlevel]" + puts "Registered sim types:" + dict for {name spec} $::sims::SIM_TYPES { + puts " $name -> [dict get $spec lang] (src/economy/sims/[dict get $spec dir]/)" + } + + puts "\nBoundedPrediction self-test:" + set bp [::sims::bounded_prediction 7.2 5.8 8.9 7.30 "4h" "amm_liquidity"] + puts " value: [dict get $bp value]" + puts " bounds: \[[dict get $bp lower_bound], [dict get $bp upper_bound]\]" + puts " confidence: [::sims::format_confidence [dict get $bp confidence]]" + puts " horizon: [dict get $bp time_horizon]" + puts " gain: [::sims::format_gain 5.8 7.2 8.9 "4h"]" + + puts "\nInvariant checks:" + if {[catch {::sims::bounded_prediction 5.0 6.0 8.0 7.0 "1h" "test"} err]} { + puts " L2 bounds check: PASS (rejected lower > value)" + } + if {[catch {::sims::bounded_prediction 5.0 4.0 8.0 11.0 "1h" "test"} err]} { + puts " Confidence range check: PASS (rejected 11.0 > 10.0)" + } + + puts "\nStatus check (no sims running):" + set st [::sims::sim_status "statistical"] + puts " statistical running: [dict get $st running]" + + puts "\nHub ready. Sims launch on M2 data feed connection." +} diff --git a/core/src/economy/sims/sociological/main.pl b/core/src/economy/sims/sociological/main.pl new file mode 100644 index 0000000..ef0630d --- /dev/null +++ b/core/src/economy/sims/sociological/main.pl @@ -0,0 +1,222 @@ +#!/usr/bin/env swipl +% +% M3b -- Sociological & population dynamics sims +% Language: Prolog (game-theoretic equilibria as constraint satisfaction) +% Protocol: line-delimited JSON on stdin/stdout to hub.tcl + +:- use_module(library(lists)). +:- use_module(library(apply)). + +% -- Behavioral archetypes ----------------------------------------------- +% Each archetype has a strategy and population share + +archetype(herd_follower). +archetype(contrarian_whale). +archetype(mev_searcher). +archetype(passive_lp). +archetype(manipulator). + +% -- Replicator dynamics -------------------------------------------------- +% dx_i/dt = x_i * [f_i(x) - phi(x)] +% Discretized: x_i(t+1) = x_i(t) + dt * x_i(t) * [f_i(x) - phi(x)] + +% Payoff matrix: row strategy vs column strategy +% Returns payoff for Row when playing against Col +payoff(herd_follower, herd_follower, 0.02). +payoff(herd_follower, contrarian_whale, -0.01). +payoff(herd_follower, mev_searcher, -0.03). +payoff(herd_follower, passive_lp, 0.01). +payoff(herd_follower, manipulator, -0.05). +payoff(contrarian_whale, herd_follower, 0.04). +payoff(contrarian_whale, contrarian_whale, -0.02). +payoff(contrarian_whale, mev_searcher, 0.01). +payoff(contrarian_whale, passive_lp, 0.02). +payoff(contrarian_whale, manipulator, -0.01). +payoff(mev_searcher, herd_follower, 0.06). +payoff(mev_searcher, contrarian_whale, 0.01). +payoff(mev_searcher, mev_searcher, -0.04). +payoff(mev_searcher, passive_lp, 0.05). +payoff(mev_searcher, manipulator, 0.02). +payoff(passive_lp, herd_follower, 0.03). +payoff(passive_lp, contrarian_whale, 0.01). +payoff(passive_lp, mev_searcher, -0.02). +payoff(passive_lp, passive_lp, 0.02). +payoff(passive_lp, manipulator, -0.04). +payoff(manipulator, herd_follower, 0.08). +payoff(manipulator, contrarian_whale, -0.03). +payoff(manipulator, mev_searcher, -0.02). +payoff(manipulator, passive_lp, 0.06). +payoff(manipulator, manipulator, -0.06). + +% Fitness of strategy I given population state Pop = [(Archetype, Share), ...] +fitness(I, Pop, F) :- + findall(Pij, ( + member((J, Xj), Pop), + payoff(I, J, Pij0), + Pij is Pij0 * Xj + ), Payoffs), + sumlist(Payoffs, F). + +% Average fitness across population +avg_fitness(Pop, Phi) :- + findall(XiFi, ( + member((I, Xi), Pop), + fitness(I, Pop, Fi), + XiFi is Xi * Fi + ), Products), + sumlist(Products, Phi). + +% One replicator step: x_i(t+dt) = x_i + dt * x_i * (f_i - phi) +replicator_step(Pop, Dt, NewPop) :- + avg_fitness(Pop, Phi), + maplist(update_share(Phi, Dt, Pop), Pop, RawPop), + normalize_pop(RawPop, NewPop). + +update_share(Phi, Dt, Pop, (I, Xi), (I, Xi1)) :- + fitness(I, Pop, Fi), + Xi1 is max(0, Xi + Dt * Xi * (Fi - Phi)). + +normalize_pop(Pop, NormPop) :- + findall(X, member((_, X), Pop), Shares), + sumlist(Shares, Total), + (Total > 0 -> + maplist(norm_share(Total), Pop, NormPop) + ; + NormPop = Pop + ). + +norm_share(Total, (I, X), (I, Xn)) :- + Xn is X / Total. + +% Run N replicator steps +replicator_evolve(Pop, _, 0, Pop) :- !. +replicator_evolve(Pop, Dt, N, FinalPop) :- + N > 0, + replicator_step(Pop, Dt, Pop1), + N1 is N - 1, + replicator_evolve(Pop1, Dt, N1, FinalPop). + +% -- Hegselmann-Krause bounded confidence --------------------------------- +% x_i(t+1) = mean({x_j : |x_j - x_i| < epsilon}) + +hk_step(Opinions, Epsilon, NewOpinions) :- + maplist(hk_update(Opinions, Epsilon), Opinions, NewOpinions). + +hk_update(AllOpinions, Epsilon, Xi, NewXi) :- + include(within_confidence(Xi, Epsilon), AllOpinions, Neighbors), + length(Neighbors, Count), + sumlist(Neighbors, Sum), + NewXi is Sum / Count. + +within_confidence(Xi, Epsilon, Xj) :- + abs(Xj - Xi) < Epsilon. + +hk_evolve(Opinions, _, 0, Opinions) :- !. +hk_evolve(Opinions, Epsilon, N, Final) :- + N > 0, + hk_step(Opinions, Epsilon, Next), + N1 is N - 1, + hk_evolve(Next, Epsilon, N1, Final). + +% -- Nash equilibrium search (constraint satisfaction) -------------------- +% For 2-player symmetric games, find mixed strategy Nash equilibria +% via support enumeration + +% Check if a mixed strategy (list of probabilities) is a Nash eq +% for a symmetric game with payoff matrix +is_nash_2p(Strategies, PayoffMatrix, Threshold) :- + length(Strategies, N), + length(PayoffMatrix, N), + expected_payoff_vec(Strategies, PayoffMatrix, ExpPayoffs), + max_list(ExpPayoffs, MaxPayoff), + forall(( + nth0(I, Strategies, Si), + nth0(I, ExpPayoffs, Ei) + ), ( + Si =:= 0 ; abs(Ei - MaxPayoff) < Threshold + )). + +expected_payoff_vec(Strat, Matrix, Payoffs) :- + maplist(expected_payoff_row(Strat), Matrix, Payoffs). + +expected_payoff_row(Strat, Row, EP) :- + maplist(mul, Strat, Row, Products), + sumlist(Products, EP). + +mul(A, B, C) :- C is A * B. + +% -- BoundedPrediction output --------------------------------------------- + +bounded_prediction(Value, Lower, Upper, Confidence, Horizon, Pred) :- + Lower =< Value, + Value =< Upper, + Confidence >= 0.0, + Confidence =< 10.0, + get_time(Now), + Timestamp is round(Now * 1000), + Pred = pred(Value, Lower, Upper, Confidence, Horizon, statistical, Timestamp). + +format_prediction(pred(V, L, U, C, H, T, _)) :- + format(" value: ~4f [~4f, ~4f]~n", [V, L, U]), + format(" confidence: ~2f/10.00~n", [C]), + format(" horizon: ~w type: ~w~n", [H, T]). + +% -- Self-test ------------------------------------------------------------- + +default_population([ + (herd_follower, 0.30), + (contrarian_whale, 0.15), + (mev_searcher, 0.10), + (passive_lp, 0.35), + (manipulator, 0.10) +]). + +run_self_test :- + prolog_flag(version, V), + format("M3b Sociological Sim -- SWI-Prolog ~w~n~n", [V]), + + format("Replicator dynamics (100 steps, dt=0.1):~n", []), + default_population(Pop0), + format(" initial: ", []), + print_pop(Pop0), + replicator_evolve(Pop0, 0.1, 100, PopFinal), + format(" final: ", []), + print_pop(PopFinal), + nl, + + format("Hegselmann-Krause (epsilon=0.2, 20 steps):~n", []), + HKInit = [0.1, 0.2, 0.25, 0.5, 0.55, 0.8, 0.85, 0.9], + format(" initial: ~w~n", [HKInit]), + hk_evolve(HKInit, 0.2, 20, HKFinal), + format(" final: ", []), + maplist(print_float, HKFinal), nl, nl, + + format("BoundedPrediction check:~n", []), + (bounded_prediction(0.35, 0.20, 0.50, 7.80, '7d', Pred) -> + format_prediction(Pred) + ; + format(" FAILED~n", []) + ), + + format("~nInvariant checks:~n", []), + (bounded_prediction(5.0, 6.0, 8.0, 7.0, '1h', _) -> + format(" L2 bounds: FAIL~n", []) + ; + format(" L2 bounds: PASS (rejected lower > value)~n", []) + ), + (bounded_prediction(5.0, 4.0, 8.0, 11.0, '1h', _) -> + format(" Confidence range: FAIL~n", []) + ; + format(" Confidence range: PASS (rejected 11.0 > 10.0)~n", []) + ), + + format("~nAll models operational.~n", []). + +print_pop([]) :- nl. +print_pop([(Name, Share)|Rest]) :- + format("~w:~3f ", [Name, Share]), + print_pop(Rest). + +print_float(X) :- format("~3f ", [X]). + +:- initialization((run_self_test, halt)). diff --git a/core/src/economy/sims/statistical/main.R b/core/src/economy/sims/statistical/main.R new file mode 100644 index 0000000..3af3bb2 --- /dev/null +++ b/core/src/economy/sims/statistical/main.R @@ -0,0 +1,275 @@ +#!/usr/bin/env Rscript +# +# M3a — Statistical & quantitative sims +# Language: R (native stats ecosystem) +# Protocol: line-delimited JSON on stdin/stdout to hub.tcl + +# -- BoundedPrediction output --------------------------------------------- + +bounded_prediction <- function(value, lower, upper, confidence, + horizon, sim_type = "statistical") { + stopifnot(lower <= value, value <= upper) + stopifnot(confidence >= 0.0, confidence <= 10.0) + list( + value = value, + lower_bound = lower, + upper_bound = upper, + confidence = confidence, + time_horizon = horizon, + sim_type = sim_type, + timestamp = as.numeric(Sys.time()) * 1000 + ) +} + +# -- Geometric Brownian Motion (Monte Carlo) ------------------------------- + +gbm_paths <- function(S0, mu, sigma, dt, n_steps, n_paths) { + Z <- matrix(rnorm(n_steps * n_paths), nrow = n_steps) + log_returns <- (mu - 0.5 * sigma^2) * dt + sigma * sqrt(dt) * Z + S <- matrix(S0, nrow = n_steps + 1, ncol = n_paths) + for (t in seq_len(n_steps)) { + S[t + 1, ] <- S[t, ] * exp(log_returns[t, ]) + } + S +} + +gbm_predict <- function(S0, mu, sigma, horizon_hours, n_paths = 10000L) { + dt <- 1 / (252 * 24) + n_steps <- as.integer(horizon_hours) + paths <- gbm_paths(S0, mu, sigma, dt, n_steps, n_paths) + final <- paths[n_steps + 1, ] + bounded_prediction( + value = median(final), + lower = quantile(final, 0.025, names = FALSE), + upper = quantile(final, 0.975, names = FALSE), + confidence = 9.50, + horizon = paste0(horizon_hours, "h") + ) +} + +# -- Heston stochastic volatility ----------------------------------------- +# dS = mu*S*dt + sqrt(nu)*S*dW^S +# dnu = kappa*(theta - nu)*dt + xi*sqrt(nu)*dW^nu +# corr(dW^S, dW^nu) = rho +# Euler-Maruyama with full truncation (nu >= 0) + +heston_paths <- function(S0, nu0, mu, kappa, theta, xi, rho, + dt, n_steps, n_paths) { + S <- matrix(S0, nrow = n_steps + 1, ncol = n_paths) + nu <- matrix(nu0, nrow = n_steps + 1, ncol = n_paths) + + for (t in seq_len(n_steps)) { + Z1 <- rnorm(n_paths) + Z2 <- rho * Z1 + sqrt(1 - rho^2) * rnorm(n_paths) + + nu_pos <- pmax(nu[t, ], 0) + sqrt_nu <- sqrt(nu_pos) + + nu[t + 1, ] <- pmax( + nu[t, ] + kappa * (theta - nu_pos) * dt + xi * sqrt_nu * sqrt(dt) * Z1, + 0 + ) + S[t + 1, ] <- S[t, ] * exp( + (mu - 0.5 * nu_pos) * dt + sqrt_nu * sqrt(dt) * Z2 + ) + } + list(S = S, nu = nu) +} + +heston_predict <- function(S0, nu0, mu, kappa, theta, xi, rho, + horizon_hours, n_paths = 5000L) { + dt <- 1 / (252 * 24) + n_steps <- as.integer(horizon_hours) + result <- heston_paths(S0, nu0, mu, kappa, theta, xi, rho, + dt, n_steps, n_paths) + final_S <- result$S[n_steps + 1, ] + final_nu <- result$nu[n_steps + 1, ] + + price_pred <- bounded_prediction( + value = median(final_S), + lower = quantile(final_S, 0.025, names = FALSE), + upper = quantile(final_S, 0.975, names = FALSE), + confidence = 9.00, + horizon = paste0(horizon_hours, "h") + ) + vol_pred <- bounded_prediction( + value = median(sqrt(final_nu)), + lower = quantile(sqrt(pmax(final_nu, 0)), 0.025, names = FALSE), + upper = quantile(sqrt(pmax(final_nu, 0)), 0.975, names = FALSE), + confidence = 8.50, + horizon = paste0(horizon_hours, "h") + ) + list(price = price_pred, volatility = vol_pred) +} + +# -- Merton jump-diffusion ------------------------------------------------ +# dS = (mu - lambda*k)*S*dt + sigma*S*dW + S*dJ +# J ~ Poisson(lambda*dt), jump size ~ LogNormal(mu_j, sigma_j) + +merton_paths <- function(S0, mu, sigma, lambda, mu_j, sigma_j, + dt, n_steps, n_paths) { + S <- matrix(S0, nrow = n_steps + 1, ncol = n_paths) + k <- exp(mu_j + 0.5 * sigma_j^2) - 1 + + for (t in seq_len(n_steps)) { + Z <- rnorm(n_paths) + N_jumps <- rpois(n_paths, lambda * dt) + J <- ifelse(N_jumps > 0, + exp(rnorm(n_paths, mu_j * N_jumps, sigma_j * sqrt(N_jumps))), + 1) + S[t + 1, ] <- S[t, ] * exp( + (mu - lambda * k - 0.5 * sigma^2) * dt + sigma * sqrt(dt) * Z + ) * J + } + S +} + +# -- HMM regime detection (3-state: Bull, Neutral, Bear) ------------------ +# r_t | s_t ~ N(mu_{s_t}, sigma_{s_t}^2) +# Forward algorithm for online filtering + +hmm_states <- c("bull", "neutral", "bear") + +hmm_filter <- function(returns, mu_vec, sigma_vec, trans_mat, init_prob) { + n <- length(returns) + K <- length(mu_vec) + alpha <- matrix(0, nrow = n, ncol = K) + + emit <- dnorm(returns[1], mu_vec, sigma_vec) + alpha[1, ] <- init_prob * emit + alpha[1, ] <- alpha[1, ] / sum(alpha[1, ]) + + for (t in 2:n) { + emit <- dnorm(returns[t], mu_vec, sigma_vec) + for (j in seq_len(K)) { + alpha[t, j] <- emit[j] * sum(alpha[t - 1, ] * trans_mat[, j]) + } + alpha[t, ] <- alpha[t, ] / sum(alpha[t, ]) + } + alpha +} + +hmm_current_regime <- function(returns, mu_vec, sigma_vec, trans_mat, + init_prob) { + alpha <- hmm_filter(returns, mu_vec, sigma_vec, trans_mat, init_prob) + last_row <- alpha[nrow(alpha), ] + state_idx <- which.max(last_row) + bounded_prediction( + value = state_idx, + lower = state_idx, + upper = state_idx, + confidence = round(max(last_row) * 10, 2), + horizon = "current" + ) +} + +# -- GARCH(1,1) volatility ------------------------------------------------ + +garch11_fit <- function(returns) { + n <- length(returns) + omega <- var(returns) * 0.05 + alpha <- 0.10 + beta <- 0.85 + sigma2 <- numeric(n) + sigma2[1] <- var(returns) + + for (t in 2:n) { + sigma2[t] <- omega + alpha * returns[t - 1]^2 + beta * sigma2[t - 1] + } + list(sigma2 = sigma2, omega = omega, alpha = alpha, beta = beta) +} + +garch_predict <- function(returns, horizon_hours) { + fit <- garch11_fit(returns) + last_var <- tail(fit$sigma2, 1) + unconditional_var <- fit$omega / (1 - fit$alpha - fit$beta) + forecast_var <- numeric(horizon_hours) + forecast_var[1] <- last_var + + for (h in 2:horizon_hours) { + forecast_var[h] <- fit$omega + + (fit$alpha + fit$beta) * forecast_var[h - 1] + } + + vol <- sqrt(mean(forecast_var)) + bounded_prediction( + value = vol, + lower = vol * 0.7, + upper = vol * 1.4, + confidence = 8.50, + horizon = paste0(horizon_hours, "h") + ) +} + +# -- Sim state (mutable, updated by ticks and calibration) ----------------- + +sim_state <- new.env(parent = emptyenv()) +sim_state$params <- list( + S0 = 1800, + mu = 0.05, + sigma = 0.60, + nu0 = 0.36, + kappa = 2.0, + theta = 0.36, + xi = 0.50, + rho = -0.70, + lambda = 5.0, + mu_j = -0.02, + sigma_j = 0.05, + hmm_mu = c(0.001, 0.0, -0.001), + hmm_sigma = c(0.01, 0.015, 0.025), + hmm_trans = matrix(c( + 0.95, 0.03, 0.02, + 0.05, 0.90, 0.05, + 0.02, 0.03, 0.95 + ), nrow = 3, byrow = TRUE), + hmm_init = c(1/3, 1/3, 1/3) +) +sim_state$price_history <- numeric(0) +sim_state$tick_count <- 0L + +# -- Self-test when run directly ------------------------------------------- + +if (!interactive() && identical(commandArgs(trailingOnly = TRUE), character(0))) { + cat("M3a Statistical Sim — R", paste(R.version$major, R.version$minor, sep = "."), "\n") + + cat("\nGBM Monte Carlo (24h, 1000 paths):\n") + bp <- gbm_predict(1800, 0.05, 0.60, 24, 1000L) + cat(sprintf(" price: %.2f [%.2f, %.2f] confidence: %.2f/10.00\n", + bp$value, bp$lower_bound, bp$upper_bound, bp$confidence)) + + cat("\nHeston SV (24h, 1000 paths):\n") + hp <- heston_predict(1800, 0.36, 0.05, 2.0, 0.36, 0.50, -0.70, 24, 1000L) + cat(sprintf(" price: %.2f [%.2f, %.2f]\n", + hp$price$value, hp$price$lower_bound, hp$price$upper_bound)) + cat(sprintf(" vol: %.4f [%.4f, %.4f]\n", + hp$volatility$value, hp$volatility$lower_bound, + hp$volatility$upper_bound)) + + cat("\nHMM regime (synthetic returns):\n") + set.seed(42) + returns <- c(rnorm(50, 0.001, 0.01), rnorm(50, -0.001, 0.025)) + rp <- hmm_current_regime(returns, + c(0.001, 0.0, -0.001), + c(0.01, 0.015, 0.025), + matrix(c(0.95,0.03,0.02, + 0.05,0.90,0.05, + 0.02,0.03,0.95), 3, byrow = TRUE), + c(1/3, 1/3, 1/3)) + cat(sprintf(" regime: %s (confidence: %.2f/10.00)\n", + hmm_states[rp$value], rp$confidence)) + + cat("\nGARCH(1,1) vol forecast:\n") + gp <- garch_predict(returns, 24) + cat(sprintf(" vol: %.6f [%.6f, %.6f]\n", + gp$value, gp$lower_bound, gp$upper_bound)) + + cat("\nMerton jump-diffusion (24h, 1000 paths):\n") + set.seed(42) + mp <- merton_paths(1800, 0.05, 0.60, 5.0, -0.02, 0.05, 1/(252*24), 24, 1000L) + final <- mp[25, ] + cat(sprintf(" median: %.2f [%.2f, %.2f]\n", + median(final), quantile(final, 0.025), quantile(final, 0.975))) + + cat("\nAll models operational.\n") +}