Thank you very much for your explanation! My code has finally debugged successfully!
I have another question: the sampling process takes a very long time. I hand-wrote an MCMC sampling procedure in Python, and running my self-written program only takes a few minutes. However, I believe using the PyMC package would definitely be much faster than the program I wrote myself. So, I’d like to ask if there’s any way to significantly increase the sampling speed?
My code is as follows:
import numpy as np
import pandas as pd
import pymc as pm
from pymc.pytensorf import collect_default_updates
import xarray as xr
import arviz as az
import pytensor
import pytensor.tensor as pt
import multiprocessing
BURNIN = 5000
N = 15000
dt = 1/240
stock_price_data = pd.read_csv("convvalue_data//2023-05-12&110043.SH.csv")
Y = np.diff(np.log(stock_price_data["convvalue"]))
T = len(Y)
V = pt.zeros(T)
def heston_dist(epsilons, V0, kappa, theta, mu, size):
def V_step(*args):
epsilons, prev_V, kappa, theta = args
new_V = epsilons[1] * pt.sqrt(prev_V * dt) + (1 - kappa * dt) * prev_V + kappa * theta * dt
return new_V, collect_default_updates(inputs=args, outputs=[new_V])
V, updates = pytensor.scan(fn=V_step, sequences=[epsilons], outputs_info=[V0], non_sequences=[kappa, theta])
Y = epsilons[:, 0] * pt.sqrt(V[:-1] * dt) - 0.5 * V[:-1] * dt + mu * dt
return Y
def calc_cov_matrix(phi, omega):
cov_matrix = pm.math.stack([
[1., phi],
[phi, phi**2 + omega]
])
return cov_matrix
with pm.Model() as mod:
# start
phi_start = np.random.normal(0, 0.01)
while True:
V0_start = np.random.normal(0.0225, 0.005)
if V0_start > 0:
break
start = {'mu': 0.1, 'kappa': 5, 'theta': 0.0225, 'phi': phi_start, 'omega': 0.02, 'V0': V0_start}
# prior distribution
mu = pm.Normal('mu', 0, sigma=1)
kappa = pm.Normal('kappa', 0, sigma=1)
theta = pm.Normal('theta', 0, sigma=1)
omega = pm.InverseGamma('omega', alpha=2, beta=0.005)
std_dev = pm.Deterministic('std_dev', omega / 2)
phi = pm.Normal('phi', 0, sigma=std_dev)
V0 = pm.TruncatedNormal('V0', mu=0.0225, sigma=0.005, lower=0)
cov_matrix = pm.Deterministic('cov_matrix', calc_cov_matrix(phi, omega))
epsilons = pm.MvNormal('epsilons', mu=np.zeros(2), cov=cov_matrix, shape=(T, 2))
y_obs = pm.MutableData('y_obs', Y)
y = pm.CustomDist('obs', epsilons, V0, kappa, theta, mu, dist=heston_dist, observed=y_obs)
idata = pm.sample_smc(draws=BURNIN+N, start=start)
Here are some implemented information:
- I’m using the Windows operating system.
- I installed PyMC using the command
pip install pymc
- The version information is as follows:
pymc.__version__ = '5.6.1'
pytensor.__version__ == '2.12.3'
