Archive Meta Data

Author(s)

Predrag Pilipovic, Adeline Samson and Susanne Ditlevsen

Title

Strang Splitting for Parametric Inference in Second-order Stochastic Differential Equations

Description

This repository is supplementary material for the paper "Strang Splitting for Parametric Inference in Second-order Stochastic Differential Equations" by Predrag Pilipovic, Adeline Samson and Susanne Ditlevsen. It provides all data generation, estimation, figures, and analysis for the Greenland Ice-Core CO2 data presented in the paper.

## Overview

This project focuses on statistical estimation for second-order stochastic differential equations (SDEs) using splitting schemes, particularly for modeling climate transitions based on ice core data. Our approach leverages the Strang splitting (SS) method to approximate the transition density between consecutive observations, which forms the basis for the pseudo-likelihood function. The exact likelihood function is often intractable, making pseudo-likelihood estimation a practical alternative. We employ maximum likelihood estimation (MLE) to estimate parameters, addressing challenges posed by hypoellipticity and partial observations.

### Model and Estimators

We analyze the Kramers oscillator, a second-order SDE used to model stochastic transitions between metastable states. The estimators compared in this study include:

- **Strang Splitting (SS) estimator**: Proposed in our work for improved accuracy and speed.
- **Euler-Maruyama (EM) estimator**: A commonly used method in applications.
- **Local Gaussian (LG) estimator** Improved EM estimator by adding an additional order of approximation in the mean of the smooth comportment and the smallest necessary orders in each of the four sub-components of the covariance matrix.
- **Local Linearization (LL) estimator**: Known for its strong performance in elliptic cases.

The true parameters of the Kramers oscillator are set to:
- \(\eta_0 = 6.5\)
- \(a_0 = 1\)
- \(b_0 = 0.6\)
- \(\sigma_0^2 = 0.1\)

### Simulation and Data

The simulation study follows these steps:

- **Trajectory simulation**: Using the EM discretization with a fine step size of \(h^{\mathrm{sim}} = 0.0001\) to ensure accuracy.
- **Sub-sampling**: Using a coarser time step \(h = 0.1\) to mitigate discretization errors.
- **Data generation**: Each trajectory consists of \(N = 5000\) data points, and we repeat simulations to obtain 500 datasets.

The study applies these estimators to proxy climate data from Greenland ice cores. Specifically, we use **calcium ion concentration (Ca²⁺) time series**, which serves as a high-resolution proxy for past climate conditions. These data originate from the Greenland Ice Core Project (GRIP) and have been processed into a uniformly spaced series using 20-year binning and averaging.

### Optimization Algorithm

For parameter estimation, we use the **resilient backpropagation algorithm (optim_rprop)** implemented in the `torch` package in R. Key optimization settings include:

- **Maximum iterations**: 2000
- **Convergence criterion**: A precision of \(10^{-5}\) for consecutive parameter estimates
- **Initial parameter values**: (-0.1, -0.1, 0.1, 0.1)

## Prerequisites

Ensure the following R packages are installed before running the project:

- `torch`: Enables automatic differentiation for fast and reliable estimators.
- `dplyr`, `tidyr`: Assists in data manipulation.
- `tictoc`: Measures code execution time.
- `reshape2`: Helps in transforming data.
- `ggplot2`, `ggpattern`, `gridExtra`, `forcats`, `ggh4x`, `ggthemes`: Enhances data visualization.
- `zoo`, `purrr`, `tidyverse`, `readxl`, `ggpubr`, `mvtnorm`: Various utilities for data handling and statistical operations.

```r
install.packages(c("torch", "dplyr", "tidyr", "tictoc", "reshape2", "ggplot2", "ggpattern",
                   "gridExtra", "forcats", "ggh4x", "ggthemes", "zoo", "purrr", "tidyverse",
                   "readxl", "ggpubr", "mvtnorm"))
```

## Usage

Run `2nd_order_main.R` to execute the whole pipeline. The dependencies are included within:

1. `2nd_order_mean_cov_functions.R`: Defines mean and covariance functions.
2. `2nd_order_objective_functions.R`: Computes objective functions for parameter estimation.
3. `2nd_order_estimator.R`: Implements optimization methods.
4. `2nd_order_data_simulator.R`: Simulates data for testing.
5. `2nd_order_results.R`: Runs estimators and saves results.
6. `2nd_order_plots.R`: Generates various plots.
7. `2nd_order_Ca_data.R`: Data analysis and estimation for Ice-Core data.
8. `2nd_order_Kramers.R`: Generates a trajectory plot for the Kramers oscillator.

Set the working directory to the project folder before execution.

### Data Regeneration and Parameter Reestimation

To regenerate data and reestimate parameters:

```r
# 1. Delete contents of '2nd_order_data' and '2nd_order_estimators' folders.
# 2. Run the required scripts in the specified order from the Usage section.
```

### Folder Structure

- `2nd_order_data`: Stores simulated data.
- `2nd_order_estimators`: Stores estimated parameters.
- `2nd_order_figures`: Stores generated figures.

### License

This project is licensed under the Creative Commons Attribution License (CC-BY).


#### Acknowledgments

This work has received funding from the European Union's Horizon 2020 research and innovation program under the Marie Skłodowska-Curie grant agreement No 956107, "Economic Policy in Complex Environments (EPOC)"; and Novo Nordisk Foundation NNF20OC0062958.

#### Contact

For any inquiries or suggestions, feel free to contact predrag@math.ku.dk