human
I am working on the file "StateSpace/Bayesian state space estimation in Python via Metropolis-Hastings.ipynb".
The first task is:
It turns out it will be convenient to write the model in terms of the precision of εε, defined to be h−1≡σ2εh−1≡σε2, and the ratio of the variances: q≡σ2η/σ2εq≡ση2/σε2 so that qh−1=σ2ηqh−1=ση2.
Then our error terms can be written:
εt∼N(0,h−1)ηt∼N(0,qh−1)
εt∼N(0,h−1)ηt∼N(0,qh−1)
And the true values are:
h−10=1/3=0.33h0−1=1/3=0.33
q=10/3=3.33q=10/3=3.33
To take a Bayesian approach to this problem, we assume that ψψ is a random variable, and we want to learn about the values of ψψ based on the data YTYT.
in fact we want a density p(ψ|YT)p(ψ|YT). To do this, we use Bayes rule to write:
p(ψ|YT)=p(YT|ψ)p(ψ)p(YT)
p(ψ|YT)=p(YT|ψ)p(ψ)p(YT)
or
p(ψ|YT)⏟posterior∝p(YT|ψ)⏟likelihoodp(ψ)⏟prior
p(ψ|YT)⏟posterior∝p(YT|ψ)⏟likelihoodp(ψ)⏟prior
The object of interest is the posterior.
to achieve it we need to specify a prior density for the unknown parameters and the likelihood function of the model.
Prior
We will use the following priors:
Precision
Since the precision must be positive, but has no theoretical upper bound, we use a Gamma prior:
h∼Gamma(αh,βh)
h∼Gamma(αh,βh)
to be specific, the density is written:
p(h)=βαhhΓ(α)hαh−1e−βhh
p(h)=βhαhΓ(α)hαh−1e−βhh
and we set the hyperparameters as αh=2,βh=2αh=2,βh=2. In this case, we have E(h)=αh/βh=1E(h)=αh/βh=1 and also E(h−1)=E(σ2ε)=1E(h−1)=E(σε2)=1.
Ratio of variances
Similarly, the ratio of variances must be positive, but has no theoretical upper bound, so we again use an (independent) Gamma prior:
q∼Gamma(αq,βq)
q∼Gamma(αq,βq)
and we set the same hyperparameters, so αq=2,βq=2αq=2,βq=2. Since E(q)=1E(q)=1, our prior is of equal variances. We then have E(σ2η)=E(qh−1)=E(q)E(h−1)=1E(ση2)=E(qh−1)=E(q)E(h−1)=1.
Initial state prior
As noted above, the Kalman filter must be initialized with μ0∼N(m0,P0)μ0∼N(m0,P0). We will use the following approximately diffuse prior:
μ0∼N(0,106)
μ0∼N(0,106)
Likelihood
For given parameters, likelihood of this model can be calculated via prediction error decomposition using an application of the Kalman filter iterations.
Posterior Simulation:...
Can you write Python code for it?
 ... Give me a the Python code you requested?
assistant
#Here is the Python code you requested.

# Priors
precision = mc.Gamma('precision', 2, 4)
ratio = mc.Gamma('ratio', 2, 1)

# Likelihood calculated using the state-space model
class LocalLevel(sm.tsa.statespace.MLEModel):
    def __init__(self, endog):
        # Initialize the state space model
        super(LocalLevel, self).__init__(endog, k_states=1,
                                         initialization='approximate_diffuse',
                                         loglikelihood_burn=1)
        
        # Initialize known components of the state space matrices
        self.ssm['design', :] = 1
        self.ssm['transition', :] = 1
        self.ssm['selection', :] = 1
    
    @property
    def start_params(self):
        return [1. / np.var(self.endog), 1.]
    
    @property
    def param_names(self):
        return ['h_inv', 'q']
    
    def update(self, params, transformed=True):
        params = super(LocalLevel, self).update(params, transformed)
        
        h, q = params
        sigma2_eps = 1. / h
        sigma2_eta = q * sigma2_eps
        
        self.ssm['obs_cov', 0, 0] = sigma2_eps
        self.ssm['state_cov', 0, 0] = sigma2_eta

# Instantiate the local level model with our simulated data
ll_mod = LocalLevel(y)

# Create the stochastic (observed) component
@mc.stochastic(dtype=LocalLevel, observed=True)
def local_level(value=ll_mod, h=precision, q=ratio):
    return value.loglike([h, q], transformed=True)

# Create the PyMC model
ll_mc = mc.Model((precision, ratio, local_level))

# Create a PyMC sample
ll_sampler = mc.MCMC(ll_mc)