<SYSTEM_TASK:>
Given the following text description, write Python code to implement the functionality described below step by step
<END_TASK>
<USER_TASK:>
Description:
 Interface to statsmodels
Step1:  Fitting a local level state space model
We assume that the observed data $y(t)$ follows
$$y(t)= f(t; \theta) + \epsilon(t),$$
where $f(t; \theta)$ is the logistic model solution. The errors are assumed to follow a local level model
Step2:  We first generate some local level noise and overlay on the ODE. Note that this noise is non-stationary, so the observed data can drift away from the actual solution for long periods.
Step3:  Note, that when wrapping the various models in this package into Pints, a useful function is the model.param_names function that provides a description of these parameters. We illustrate this below
Step4:  Wrapping the a state space model in a Pints log-likelihood. Here, ss_model is a string corresponding to the list of prebuilt state space models here; ss_n_params is the number of parameters that the particular time series model involves.
Step5:  Instantiate the model, log-posterior etc., then run MCMC.
Step6:  Show traces.
Step7:  Look at summary stats
<ASSISTANT_TASK:>
Python Code:
import pints
import pints.toy as toy
import pints.plot
import numpy as np
import matplotlib.pyplot as plt
Explanation: Interface to statsmodels: state space time series models
This notebook provides a short exposition of how it is possible to interface with the cornucopia of time series models provided by the statsmodels package. In this notebook, we illustrate how to fit the logistic ODE model, where the errors are described by state space models.
End of explanation
from statsmodels.tsa.statespace.structural import UnobservedComponents
Explanation: Fitting a local level state space model
We assume that the observed data $y(t)$ follows
$$y(t)= f(t; \theta) + \epsilon(t),$$
where $f(t; \theta)$ is the logistic model solution. The errors are assumed to follow a local level model:
$$
\begin{align}
&\epsilon(t) = \mu(t) + \nu(t),\
&\mu(t) = \mu(t-1) + \eta(t),
\end{align}
$$
as described here (also see this link for a description of the huge family of state space models available in this package). Here, $\nu(t) \sim \mathcal{N}(0, \sigma_\nu)$ and $\eta(t) \sim \mathcal{N}(0, \sigma_\eta)$.
The two parameters of the error process are the variances: $\sigma_\nu^2$ and $\sigma_\eta^2$.
End of explanation
import scipy.stats
# Load a forward model
model = toy.LogisticModel()
# Create some toy data
real_parameters = [0.015, 500]
times = np.linspace(0, 1000, 1000)
org_values = model.simulate(real_parameters, times)
# function to generate local level errors
def local_level(n, sigma_nu, sigma_eta):
    nu = scipy.stats.norm.rvs(0, sigma_nu, n)
    eta = scipy.stats.norm.rvs(0, sigma_eta, n)
    epsilon = np.zeros(n)
    mu = np.zeros(n)
    mu[0] = eta[0]
    for i in range(1, n):
        mu[i] = mu[i - 1] + eta[i]
    for i in range(n):
        epsilon[i] = mu[i] + nu[i]
    return epsilon
sigma_nu = 10
sigma_eta = 1
errors = local_level(len(org_values), sigma_nu, sigma_eta)
values = org_values + errors
# Show the noisy data
plt.figure()
plt.plot(times, org_values)
plt.plot(times, values)
plt.xlabel('time')
plt.ylabel('y')
plt.legend(['true', 'observed'])
plt.show()
Explanation: We first generate some local level noise and overlay on the ODE. Note that this noise is non-stationary, so the observed data can drift away from the actual solution for long periods.
End of explanation
sol = model.simulate([0.15, 500], times)
model = UnobservedComponents(endog=values, level='llevel', exog=sol)
model.param_names
Explanation: Note, that when wrapping the various models in this package into Pints, a useful function is the model.param_names function that provides a description of these parameters. We illustrate this below: note that beta.x1 is the coefficient on the ODE solution, so needs to be supplied as '1' in the wrapper we show later.
End of explanation
class StateSpaceLogLikelihood(pints.ProblemLogLikelihood):
    def __init__(self, problem, ss_model, ss_n_params):
        super(StateSpaceLogLikelihood, self).__init__(problem)
        self._nt = len(self._times) - 1
        self._no = problem.n_outputs()
        self._ss_model = ss_model
        self._ss_n_params = ss_n_params
        self._n_parameters = problem.n_parameters() + ss_n_params * self._no
        self._m = self._ss_n_params * self._no
        
    def __call__(self, x):
        # converting to list makes it easier to append
        # nuisance params
        x = x.tolist()
        m = self._m
        # extract noise model params
        parameters = x[-m:]
        sol = self._problem.evaluate(x[:-m])
        model = UnobservedComponents(endog=problem.values(),
                                     level=self._ss_model, exog=sol)
        
        # add nuisance parameter at end: the coefficient
        # on the ODE model solution (which should be 1)
        return model.loglike(parameters + [1])
Explanation: Wrapping the a state space model in a Pints log-likelihood. Here, ss_model is a string corresponding to the list of prebuilt state space models here; ss_n_params is the number of parameters that the particular time series model involves.
End of explanation
model = toy.LogisticModel()
# Create an object with links to the model and time series
problem = pints.SingleOutputProblem(model, times, values)
# Create a log-likelihood function
log_likelihood = StateSpaceLogLikelihood(problem, 'llevel', 2)
# Create a uniform prior over both the parameters and the new noise variable
log_prior = pints.UniformLogPrior(
    [0.00, 400, 0, 0],
    [0.05, 600, 3000, 3000],
)
# Create a posterior log-likelihood (log(likelihood * prior))
log_posterior = pints.LogPosterior(log_likelihood, log_prior)
# Choose starting points for 3 mcmc chains
real_parameters = np.array([0.015, 500] + [sigma_nu**2, sigma_eta**2])
xs = [
    real_parameters * 1.05,
    real_parameters * 1,
    real_parameters * 1.025
]
# Create mcmc routine
mcmc = pints.MCMCController(log_posterior, 3, xs, method=pints.HaarioBardenetACMC)
# Add stopping criterion
mcmc.set_max_iterations(4000)
# Disable logging
mcmc.set_log_to_screen(False)
# Run!
print('Running...')
chains = mcmc.run()
print('Done!')
Explanation: Instantiate the model, log-posterior etc., then run MCMC.
End of explanation
pints.plot.trace(chains,
                 parameter_names=[r'$r$', r'$k$', r'$\sigma^2_{\nu}$', r'$\sigma^2_{\eta}}$'],
                 ref_parameters=real_parameters)
plt.show()
Explanation: Show traces.
End of explanation
results = pints.MCMCSummary(chains=chains, 
                            parameter_names=["r", "k", "sigma^2_nu", "sigma^2_eta"])
print(results)
Explanation: Look at summary stats: all parameters look similar to their true values.
End of explanation
<END_TASK>