SPX Parameter Estimation & Convergence

Performance, Convergence, and Likelihood Validation: R pomp vs. pypomp

Published

September 2, 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 parameter convergence and likelihood evaluation of the S&P 500 (SPX) stochastic volatility model used in Sun (2024) using the pypomp package (on CPU and GPU) and the R pomp package (on CPU).

The SPX model is a stochastic volatility model with the following parameter set:

  • \(\mu\): Expected rate of return of the asset.
  • \(\kappa\): Mean rate of reversion in variance.
  • \(\theta\): Long term average variance.
  • \(\xi\): Volatility coefficient of volatility.
  • \(\rho\): Leverage effect parameter (correlation between return and volatility shocks).
  • \(V_0\): Initial variance.

This benchmark compares runtimes, convergence trajectories, and likelihood distributions to verify the mathematical correctness and efficiency of pypomp.

Benchmark Settings & Environment

The table below summarizes algorithmic parameters, software environments, and compute hardware recorded in latest.json for each configuration.

Show Code 
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)
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
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>"))
Show Code 
import json
import gc
import numpy as np
import pandas as pd
from plotnine import (
    ggplot, aes, geom_line, geom_ribbon, geom_density, facet_wrap,
    labs, theme, guides, guide_legend
)

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

import pypomp as pp
import jax
Show Code 
GPU_JSON_PATH = os.path.join("results", "gpu", "latest.json")
CPU_JSON_PATH = os.path.join("results", "cpu", "latest.json")

GPU_CSV_PATH = os.path.join("results", "gpu", "results.csv")
CPU_CSV_PATH = os.path.join("results", "cpu", "results.csv")

GPU_TRACES_PATH = os.path.join("results", "gpu", "traces.csv.gz")
CPU_TRACES_PATH = os.path.join("results", "cpu", "traces.csv.gz")

gpu_json = json.load(open(GPU_JSON_PATH)) if os.path.exists(GPU_JSON_PATH) else {}
cpu_json = json.load(open(CPU_JSON_PATH)) if os.path.exists(CPU_JSON_PATH) else {}

LL_frame_gpu, traces_gpu = ru.load_results_and_traces(GPU_CSV_PATH, GPU_TRACES_PATH)
LL_frame_cpu, traces_cpu = ru.load_results_and_traces(CPU_CSV_PATH, CPU_TRACES_PATH)

Parameter Convergence & Traces (IF2)

Comparison of Python and R Parameter Traces

Show Code 
last_iter = traces_gpu["iteration"].max() if "iteration" in traces_gpu.columns else traces_gpu["iter"].max()
gpu_pfilter_last = traces_gpu[(traces_gpu["iteration" if "iteration" in traces_gpu.columns else "iter"] == last_iter) & (traces_gpu["method"] == "pfilter")]
best_rep = gpu_pfilter_last.groupby("theta_idx" if "theta_idx" in gpu_pfilter_last.columns else "rep")["logLik"].first().idxmax()

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)")
del traces_gpu

best_gpu_pfilter = traces_gpu_t[
    (traces_gpu_t["iter"] == last_iter) & 
    (traces_gpu_t["rep"] == best_rep) & 
    (traces_gpu_t["method"] == "pfilter")
].copy()

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

traces_long_R = pd.read_csv(os.path.join(R_REF, "mif_traces.csv.gz"))
traces_long_R["iter"] = (traces_long_R["iteration"] + 1).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()

Here, we overlay the trajectories from R, Python (CPU), and Python (GPU) to compare the search paths directly.

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}))
)

Expectation

Since pypomp is designed to be a direct Python/JAX implementation of R’s pomp package, the paths taken by the individual replicates in Python (both CPU and GPU) should be similar to those in R, overlapping substantially.

Comparison of Python and R Parameter Trace Quantiles

To abstract away the noise of individual runs, this plot compares the 10th and 90th percentile bounds of the parameters across all runs.

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
)

Expectation

The 10th and 90th percentile ribbons provide a robust view of convergence. The ribbons for Python (GPU), Python (CPU), and R should overlap closely throughout the iterations.

Parameter Estimates Distribution

Comparison of Final Parameter Estimates

This plot compares the distribution of final estimated parameters across replicates.

Show Code 
max_iters = traces_pyr.groupby(["language", "quantity"], observed=True)["iter"].transform("max")
final_traces = traces_pyr[traces_pyr["iter"] == max_iters].copy()
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
)

Expectation

The density of the final parameter values across replicates should overlap heavily for each configuration.

Comparison of Starting Parameter Values (Sanity Check)

Show Code 
init_traces = traces_pyr[
    (~traces_pyr["quantity"].isin(["logLik", "loglik"])) & 
    (traces_pyr["iter"] == traces_pyr.groupby(["language", "quantity"], observed=True)["iter"].transform("min"))
]

(
    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
)

Expectation

Because CPU, GPU, and R runs are initialized using the same random starting box, the starting densities should be very similar across configurations.

Log-Likelihood Comparisons

Summary of Final Log-Likelihood Estimates

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

Python (GPU) Top 10:

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

Python (CPU) Top 10:

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 Top 10:

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

Density Plot Comparison (LL > 11,800)

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
)

Expectation

Final log-likelihood estimates from each configuration should overlap heavily.

References

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