SPX Parameter Estimation & Convergence

IF2 global search: R pomp vs. pypomp

Published

October 16, 2026

Show Code

import os
import sys
sys.path.append("..")
import report_utils as ru
from IPython.display import display, HTML

display(HTML(f"<div style='margin-bottom: 20px;'>{ru.nav_bar('estimation')}</div>"))

Estimation | Timing & Throughput | Scaling | Likelihood Evaluation

Introduction

This report compares IF2 global searches on the S&P 500 (SPX) stochastic volatility model of Sun (2024) in pypomp (GPU and CPU) and R pomp. Each implementation starts from points drawn from the same box and uses the same number of iterations and perturbation sizes. We compare the starting values, the final parameter estimates and the final log-likelihoods across starts.

The model has six parameters:

$\mu$: expected rate of return of the asset.
$\kappa$: mean rate of reversion in variance.
$\theta$: long-term average variance.
$\xi$: volatility of volatility.
$\rho$: leverage effect, the correlation between return and volatility shocks.
$V_0$: initial variance.

Parameters are compared on the transformed scale (log for $\mu$, $\kappa$, $\theta$, $\xi$ and $V_0$, and $\log\frac{1+\rho}{1-\rho}$ for $\rho$).

Show Code

import gc
import numpy as np
import pandas as pd
import checks
from plotnine import (
    ggplot, aes, geom_line, geom_ribbon, geom_density, facet_wrap,
    labs, guides, guide_legend
)

os.environ["JAX_PLATFORMS"] = "cpu"

PLATFORMS = {
    "R pomp (36 cores)": os.path.join("results", "R"),
    "pypomp (GPU)": os.path.join("results", "gpu"),
    "pypomp (CPU, 36 cores)": os.path.join("results", "cpu"),
}
runs = ru.load_timing_data(PLATFORMS)

LL_frame_gpu, traces_gpu = ru.load_results_and_traces(
    os.path.join("results", "gpu", "results.csv"), os.path.join("results", "gpu", "traces.csv.gz")
)
LL_frame_cpu, traces_cpu = ru.load_results_and_traces(
    os.path.join("results", "cpu", "results.csv"), os.path.join("results", "cpu", "traces.csv.gz")
)

R_REF = os.path.join("results", "R")
LL_df_R = pd.read_csv(os.path.join(R_REF, "pfilter_logliks.csv"))[["logLik", "se"]].rename(columns={"logLik": "est"})

traces_gpu_t = ru.process_and_transform_traces(traces_gpu, "python (GPU)")
traces_cpu_t = ru.process_and_transform_traces(traces_cpu, "python (CPU)")
del traces_gpu, traces_cpu

traces_long_R = ru.align_pomp_mif_loglik(
    pd.read_csv(os.path.join(R_REF, "mif_traces.csv.gz")), ["replicate"]
)
traces_long_R["iter"] = traces_long_R["iteration"].astype(np.int16)
traces_long_R["rep"] = traces_long_R["replicate"].astype(np.int16)
traces_long_R = traces_long_R.drop(columns=["iteration", "replicate"])
traces_R_t = ru.process_and_transform_traces(traces_long_R, "R")
del traces_long_R

traces_pyr = pd.concat([traces_gpu_t, traces_cpu_t, traces_R_t], ignore_index=True)
del traces_gpu_t, traces_cpu_t, traces_R_t
_ = gc.collect()

is_param = ~traces_pyr["quantity"].isin(["logLik", "loglik"])
max_iters = traces_pyr.groupby(["language", "quantity"], observed=True)["iter"].transform("max")
min_iters = traces_pyr.groupby(["language", "quantity"], observed=True)["iter"].transform("min")
final_traces = traces_pyr[traces_pyr["iter"] == max_iters].copy()
init_traces = traces_pyr[is_param & (traces_pyr["iter"] == min_iters)]

Show Code

missing = [(label, r["dir"]) for label, r in runs.items() if not r["available"]]
for label, path in missing:
    display(HTML(f"<div class='alert alert-warning'><strong>Missing results for {label}:</strong> Expected at <code>{path}</code>.</div>"))

Automated checks

Show Code

PY_LANGS = {"python (GPU)": "pypomp (GPU)", "python (CPU)": "pypomp (CPU)"}
PARAMS = ["mu", "kappa", "theta", "xi", "rho", "V_0"]
PY_LL = {"python (GPU)": LL_frame_gpu["LL"], "python (CPU)": LL_frame_cpu["LL"]}

# One row per start and parameter; at the last iteration pypomp also records a pfilter row.
by_start = {
    "Final estimates": final_traces[final_traces["quantity"].isin(PARAMS)].drop_duplicates(["language", "quantity", "rep"]),
    "Starting values": init_traces.drop_duplicates(["language", "quantity", "rep"]),
}


def samples(df, lang):
    return [
        (q, df[(df["language"] == lang) & (df["quantity"] == q)]["value_T"],
         df[(df["language"] == "R") & (df["quantity"] == q)]["value_T"])
        for q in PARAMS
    ]


M_TESTS = len(PY_LANGS) * (1 + 2 * len(by_start) * len(PARAMS))
rows, details = [], []
for lang, label in PY_LANGS.items():
    compared = f"{label} vs R pomp"
    rows += checks.ks_check(PY_LL[lang], LL_df_R["est"], "Final log-likelihood (KS)", compared, M_TESTS)
    for check, df in by_start.items():
        r, det = checks.group_checks(samples(df, lang), check, compared, M_TESTS)
        rows += r
        details.append(det.assign(check=check))

level = checks.run_level(runs)
rows = checks.gate_smoke(rows, level)
checks.save(rows, level, pypomp_version=runs["pypomp (GPU)"]["meta"].get("pypomp_version"))
display(HTML(checks.table_html(rows, M_TESTS)))
display(HTML(checks.details_html(pd.concat(details, ignore_index=True))))

These checks are computed from the results each time the report is rendered. They only compare summary statistics, so the plots can show problems that the checks miss. Visual review lists what the plots should show.

Check	Compared	Result	Verdict
Final log-likelihood (KS)	pypomp (GPU) vs R pomp	D = 0.0528, p = 0.698	PASS
Final estimates: means	pypomp (GPU) vs R pomp	largest \|z\| = 1.85 (xi) of 6	PASS
Final estimates: distributions (KS)	pypomp (GPU) vs R pomp	smallest p = 0.0234 (kappa) of 6	PASS
Starting values: means	pypomp (GPU) vs R pomp	largest \|z\| = 1.44 (kappa) of 6	PASS
Starting values: distributions (KS)	pypomp (GPU) vs R pomp	smallest p = 0.226 (kappa) of 6	PASS
Final log-likelihood (KS)	pypomp (CPU) vs R pomp	D = 0.0694, p = 0.351	PASS
Final estimates: means	pypomp (CPU) vs R pomp	largest \|z\| = 0.69 (theta) of 6	PASS
Final estimates: distributions (KS)	pypomp (CPU) vs R pomp	smallest p = 0.401 (kappa) of 6	PASS
Starting values: means	pypomp (CPU) vs R pomp	largest \|z\| = 1.44 (kappa) of 6	PASS
Starting values: distributions (KS)	pypomp (CPU) vs R pomp	smallest p = 0.226 (kappa) of 6	PASS

PASS when |z| < 3, CHECK when 3 ≤ |z| < 5, and FAIL when |z| ≥ 5. KS p-values use the equivalent thresholds. Both are Bonferroni-adjusted over the 50 tests that receive a verdict.

Results for each parameter

compared	name	difference	SE	z	KS D	KS p	check
pypomp (GPU) vs R pomp	mu	0.06845	0.3175	0.2156	0.03889	0.9488	Final estimates
pypomp (GPU) vs R pomp	kappa	-0.03197	0.05339	-0.5987	0.1111	0.02342	Final estimates
pypomp (GPU) vs R pomp	theta	-0.01615	0.01699	-0.9502	0.04722	0.8177	Final estimates
pypomp (GPU) vs R pomp	xi	0.3588	0.1936	1.853	0.08056	0.1934	Final estimates
pypomp (GPU) vs R pomp	rho	-0.05168	0.2042	-0.253	0.1056	0.03616	Final estimates
pypomp (GPU) vs R pomp	V_0	0.007629	0.07319	0.1042	0.09167	0.0971	Final estimates
pypomp (GPU) vs R pomp	mu	0.07016	0.06724	1.043	0.04167	0.9141	Starting values
pypomp (GPU) vs R pomp	kappa	-0.1062	0.07361	-1.442	0.07778	0.2265	Starting values
pypomp (GPU) vs R pomp	theta	0.015	0.01993	0.7527	0.05556	0.6357	Starting values
pypomp (GPU) vs R pomp	xi	-0.004986	0.08083	-0.06169	0.03889	0.9488	Starting values
pypomp (GPU) vs R pomp	rho	-0.07767	0.09135	-0.8503	0.05	0.7599	Starting values
pypomp (GPU) vs R pomp	V_0	0.07647	0.07544	1.014	0.03611	0.9733	Starting values
pypomp (CPU) vs R pomp	mu	-0.2207	0.3265	-0.6758	0.03889	0.9488	Final estimates
pypomp (CPU) vs R pomp	kappa	0.01364	0.05507	0.2477	0.06667	0.4009	Final estimates
pypomp (CPU) vs R pomp	theta	-0.01303	0.01883	-0.6919	0.06111	0.5127	Final estimates
pypomp (CPU) vs R pomp	xi	-0.05117	0.2362	-0.2166	0.03889	0.9488	Final estimates
pypomp (CPU) vs R pomp	rho	-0.075	0.1966	-0.3816	0.05377	0.6495	Final estimates
pypomp (CPU) vs R pomp	V_0	0.04399	0.07528	0.5844	0.05833	0.5733	Final estimates
pypomp (CPU) vs R pomp	mu	0.07017	0.06724	1.044	0.04167	0.9141	Starting values
pypomp (CPU) vs R pomp	kappa	-0.1062	0.07361	-1.442	0.07778	0.2265	Starting values
pypomp (CPU) vs R pomp	theta	0.015	0.01993	0.7528	0.05556	0.6357	Starting values
pypomp (CPU) vs R pomp	xi	-0.004988	0.08083	-0.06171	0.03889	0.9488	Starting values
pypomp (CPU) vs R pomp	rho	-0.07767	0.09135	-0.8503	0.05	0.7599	Starting values
pypomp (CPU) vs R pomp	V_0	0.07646	0.07544	1.014	0.03611	0.9733	Starting values

The mean test is not applied to the final log-likelihoods, since a few starts in each implementation end far below the rest and dominate the mean.

Visual review

Parameter traces: the trajectories from pypomp (GPU and CPU) and pomp should overlap substantially.
Trace quantiles: the 10th to 90th percentile ribbons should overlap closely throughout the iterations.
Final parameter estimates: the densities should overlap heavily for each parameter.
Starting values: the densities should be very similar across configurations.
Final log-likelihoods: the densities should overlap heavily.

Summary statistics

Show Code

ll_summary_table = pd.DataFrame({
    "Configuration": ["pypomp (GPU)", "pypomp (CPU)", "R pomp"],
    "Min": [LL_frame_gpu["LL"].min(), LL_frame_cpu["LL"].min(), LL_df_R["est"].min()],
    "Median": [LL_frame_gpu["LL"].median(), LL_frame_cpu["LL"].median(), LL_df_R["est"].median()],
    "Mean": [LL_frame_gpu["LL"].mean(), LL_frame_cpu["LL"].mean(), LL_df_R["est"].mean()],
    "Max": [LL_frame_gpu["LL"].max(), LL_frame_cpu["LL"].max(), LL_df_R["est"].max()],
    "SD": [LL_frame_gpu["LL"].std(), LL_frame_cpu["LL"].std(), LL_df_R["est"].std()]
})

display(HTML(ll_summary_table.round(2).to_html(classes="table table-striped table-hover", index=False)))

Configuration	Min	Median	Mean	Max	SD
pypomp (GPU)	1.082221e+04	11841.59	1.179416e+04	11852.42	2.121800e+02
pypomp (CPU)	1.082112e+04	11840.77	1.176843e+04	11851.09	2.595000e+02
R pomp	-7.886362e+28	11841.38	-2.190656e+26	11849.36	4.156478e+27

Top 10 starts for each configuration

pypomp (GPU):

Show Code

display(HTML(LL_frame_gpu.sort_values(by="LL", ascending=False).head(10).to_html(classes="table table-striped table-hover", index=False)))

LL	sd
11852.418993	2.545271
11851.742004	3.526382
11850.029043	2.271595
11849.986510	2.180859
11849.501113	0.435737
11849.383609	1.679806
11849.146096	0.632354
11848.871077	0.640241
11848.636547	1.874580
11848.632121	0.660764

pypomp (CPU):

Show Code

display(HTML(LL_frame_cpu.sort_values(by="LL", ascending=False).head(10).to_html(classes="table table-striped table-hover", index=False)))

LL	sd
11851.094203	1.690370
11849.907911	1.690701
11849.718760	1.377921
11849.168175	2.620740
11849.102231	0.581125
11848.940299	1.505824
11848.828586	0.342690
11848.768941	0.922149
11848.759355	0.435271
11848.596543	0.499511

R pomp:

Show Code

display(HTML(LL_df_R.sort_values(by="est", ascending=False).head(10).to_html(classes="table table-striped table-hover", index=False)))

est	se
11849.357528	0.361169
11849.291865	2.371300
11849.027270	0.608924
11848.946284	0.777941
11848.898417	1.179798
11848.837613	0.428583
11848.719805	0.790108
11848.647530	0.494265
11848.534922	2.326628
11848.367994	1.443287

Parameter traces

Show Code

traces_plot_df = traces_pyr[(traces_pyr["quantity"] != "logLik") & (traces_pyr["value_T"].between(-12, 12))].copy()
traces_plot_df["group_id"] = (traces_plot_df["rep"].astype(str) + "_" + traces_plot_df["language"].astype(str)).astype("category")

(
    ggplot(traces_plot_df, aes(x="iter", y="value_T", group="group_id", color="language"))
    + geom_line(alpha=0.15)
    + facet_wrap("~quantity", scales="free")
    + labs(
        title="Comparison of Python and R Parameter Traces",
        subtitle="Overlay of all replicate trajectories (transformed values)",
        x="Iteration",
        y="Transformed Value",
        color="Configuration"
    )
    + ru.scale_color_premium()
    + ru.theme_premium
    + guides(color=guide_legend(override_aes={"alpha": 1}))
)

Trace quantiles

Show Code

quant_df = (
    traces_plot_df.groupby(["iter", "quantity", "language"], observed=True)["value_T"]
    .quantile([0.1, 0.9])
    .unstack(level=-1)
    .reset_index()
    .rename(columns={0.1: "ymin", 0.9: "ymax"})
)
quant_df["ymid"] = (quant_df["ymin"] + quant_df["ymax"]) / 2

(
    ggplot(quant_df, aes(x="iter", ymin="ymin", ymax="ymax", fill="language", color="language"))
    + geom_ribbon(alpha=0.3, color="none")
    + geom_line(aes(y="ymid"), linetype="dashed", alpha=0.7)
    + facet_wrap("~quantity", scales="free")
    + labs(
        title="Comparison of Python and R Parameter Trace Quantiles",
        subtitle="Shaded regions represent the 10th to 90th percentile interval; dashed line is the midpoint",
        x="Iteration",
        y="Transformed Value Range",
        fill="Configuration",
        color="Configuration"
    )
    + ru.scale_fill_premium()
    + ru.scale_color_premium()
    + ru.theme_premium
)

Final parameter estimates

Show Code

final_traces_zoom = final_traces[final_traces["value_T"].between(-10, 10)]

(
    ggplot(final_traces_zoom, aes(x="value_T", fill="language", color="language"))
    + geom_density(alpha=0.3)
    + facet_wrap("~quantity", scales="free")
    + labs(
        title="Comparison of Final Parameter Estimates",
        subtitle="Density distribution of parameters at the final IF2 iteration (transformed, zoomed)",
        x="Transformed Parameter Value",
        y="Density",
        fill="Configuration",
        color="Configuration"
    )
    + ru.scale_fill_premium()
    + ru.scale_color_premium()
    + ru.theme_premium
)

Starting values

Show Code

(
    ggplot(init_traces, aes(x="value_T", fill="language", color="language"))
    + geom_density(alpha=0.3)
    + facet_wrap("~quantity", scales="free")
    + labs(
        title="Comparison of Initial Parameter Distributions",
        subtitle="Density of starting parameter configurations across restarts (transformed, zoomed)",
        x="Transformed Parameter Value",
        y="Density",
        fill="Configuration",
        color="Configuration"
    )
    + ru.scale_fill_premium()
    + ru.scale_color_premium()
    + ru.theme_premium
)

Final log-likelihoods

Starts with a final log-likelihood at or below 11,800 are left out of the plot.

Show Code

ll_density_df = pd.concat([
    pd.DataFrame({"LL": LL_frame_gpu["LL"], "language": "python (GPU)"}),
    pd.DataFrame({"LL": LL_frame_cpu["LL"], "language": "python (CPU)"}),
    pd.DataFrame({"LL": LL_df_R["est"], "language": "R"})
], ignore_index=True)

ll_density_filtered = ll_density_df[ll_density_df["LL"] > 11800]

(
    ggplot(ll_density_filtered, aes(x="LL", fill="language", color="language"))
    + geom_density(alpha=0.3)
    + labs(
        title="Comparison of Final Log-Likelihood Estimates",
        subtitle="Density of estimated log-likelihoods at final parameter sets (LL > 11,800)",
        x="Log-Likelihood (LL)",
        y="Density",
        fill="Configuration",
        color="Configuration"
    )
    + ru.scale_fill_premium()
    + ru.scale_color_premium()
    + ru.theme_premium
)

Appendix: settings and environment

Show Code

display(HTML(ru.build_settings_comparison_html(runs, is_panel=False)))

Setting / Parameter	R pomp (36 cores)	pypomp (GPU)	pypomp (CPU, 36 cores)
Algorithmic & Workload Settings
Run Level	4	4	4
Starting Searches ($N_{starts}$)	360	360	360
IF2 Iterations ($N_{iter}$)	200	200	200
IF2 Particles ($N_p$)	1,000	1,000	1,000
Pfilter Particles ($N_{p,eval}$)	1,000	1,000	1,000
Evaluation Replicates ($N_{reps}$)	24	24	24
Random Seed	—	631409	631409
Software & Environment
Pomp Framework	pomp 6.3	pypomp 1.0.0rc1	pypomp 1.0.0rc1
Backend / Engine	R 4.4.0	JAX 0.11.1	JAX 0.11.1
Quant Git Commit	`8853fc5`	`630cc23`	`630cc23`
Run Timestamp	2026-08-05 22:40:59	2026-09-02 18:33:59	2026-09-02 19:12:02
Hardware & Compute
Compute Device	Intel(R) Xeon(R) Gold 6154 CPU @ 3.00GHz (36 cores)	Tesla V100-PCIE-16GB (1 GPU)	Intel(R) Xeon(R) Gold 6254 CPU @ 3.10GHz (36 cores)
Slurm Partition	standard	gpu	standard
Slurm Job ID	56537485	59673765	59673766

References

Sun, Weizhe. 2024. “Model Based Inference of Stochastic Volatility via Iterated Filtering.” Undergraduate Honors Thesis, University of Michigan.