fix AR(1) log likelihood

[eph]

Fixes #2284

I've never contributed to an open source project, so I'm not sure if a pull request if overkill here.

[junpenglao]

Welcome! No contribution is too small ;-)
Another point of reference is the implementation in statsmodels: doc; code.

[twiecki]

Yeah, this looks like a bug to me.

[twiecki]

pep8 wants spaces around math operators.

[twiecki]

But yeah, this test is a bit odd as it basically repcliates the likelihood. statsmodels is an idea but we don't want that as a dependency, but could be optional. @junpenglao?

[junpenglao]

I dont think we should use statsmodels in the test. I am thinking more rewriting our AR1 logp similar to the (exact) unconditional maximum likelihood as in statsmodels. We should take this opportunity to extend our AR1 model into AR(n).

[eph]

• Sorry about the test. I didn't know the best way to write one without introducing additional dependencies.
• As it written now, the AR(1) likelihood calculation is conditional on $x_0 = 0$, which seems a bit arbitrary to me. It seems more natural to conditional on the initial observation.
• You could use the unconditional distribution to calculate the likelihood of the initial observation (as in stats models), but this only makes sense if the process is stationary ($ |k| < 1 $). Maybe something like init in GaussianRandomWalk?
• I am trying to write a more general AR RV, but I'm getting a bit tripped up by creating the matrix of lagged observations in theano.

[twiecki]

@eph Any progress?

[junpenglao]

@twiecki I think this is waiting for after 3.1rc4

[twiecki]

@junpenglao Definitely, but can start moving some things into the 3.2 branch when they are done.

[ferrine]

Is it ready?

[junpenglao]

I think we can merge this and improve the test later.

[ferrine]

Don't you think we'll forget about test. Or do we need test here?

[junpenglao]

@eph wrote a test here but it basically replicates the likelihood (see @twiecki's comment above).
But it actually wont be a big problem if we are going to rewrite the logp later closer to the form in statsmodels

[junpenglao]

pinging @eph

[eph]

I have written a more general (p lag) autoregressive model. Should I replace AR1 with this class (after tests are written of course). What do you want the tests written against? I know statsmodels was mentioned, but 1) the loglikelihood calculation there is not completely clear to me and 2) I thought we didn't want to keep that dependency.

class AR(distribution.Continuous):
    R"""
    Autoregressive process with p lags.

    .. math::

       x_t = \rho_0 + \rho_1 x_{t-1} + \ldots + \rho_p x_{t-p} + \epsilon_t,
       \epsilon_t \sim N(0,\sigma^2)

    
    The innovation can be parameterized either in terms of precision
    or standard deviation. The link between the two parametrizations is
    given by

    .. math::

       \tau = \dfrac{1}{\sigma^2}


    Parameters
    ----------
    rho : tensor
        Vector of autoregressive coefficients. 

    sd : float 
        Standard deviation of innovation (sd > 0).

    tau : float
        Precision of innovation (tau > 0).

    constant: bool (optional, default = False)
        Whether to include a constant. 
    """

    def __init__(self, rho, sd=None, tau=None,
                 constant=False, *args, **kwargs):
        
        super(AR, self).__init__(*args, **kwargs)
        tau, sd = get_tau_sd(tau=tau, sd=sd)
        self.sd = tt.as_tensor_variable(sd)
        self.tau = tt.as_tensor_variable(tau)

        self.mean = tt.as_tensor_variable(0.)

        rho = tt.as_tensor_variable(rho,ndim=1)
        if constant:
            self.p = rho.shape[0] - 1
        else:
            self.p = rho.shape[0]

        #if self.p <= 0:
        #    raise ValueError('too few lags.')

        self.constant = constant
        self.rho = rho
        

    def logp(self, value):

        y = value[self.p:]
        results,_ = scan(lambda l, obs, p: obs[p-l:-l],
                         outputs_info=None,sequences=[tt.arange(1,self.p+1)],
                         non_sequences=[value, self.p])
        x = tt.stack(results)

        if self.constant:
            y = y - self.rho[0]
            eps = y - self.rho[1:].dot(x)
        else:
            eps = y - self.rho.dot(x)

        innov_like = Normal.dist(mu=0.0, tau=self.tau).logp(eps)
        return tt.sum(innov_like)

[twiecki]

Yes, general AR would be great. Tests are difficult here, I wouldn't block on it. Perhaps just a regression test with fixed values? Or maybe a simple case where the result can be confirmed by hand?

[junpenglao]

There are conflicts @eph

[ferrine]

We usually do not use blank lines between

[junpenglao]

This should no appear.

[junpenglao]

This should not be deleted.

[junpenglao]

hi @eph, could you please edit the PR following the review? We should merge it after test pass.

[ferrine]

The notebook needs to be refined

[junpenglao]

Hi @eph, hope you dont mind I step in to fix the conflict and the notebook.

As for this PR, how should we process? The test should be further refined, and the general AR logp could be improved as well (without using theano.scan). But in my opinion, the pm.timeseries needs some refactoring and further documentation, so it might be better to address these in a separate PR.

[junpenglao]

OK, after discuss with @twiecki i am merging this now. Further work on AR is referenced in #2491.

[junpenglao]

Thanks for the contribution @eph!

[eph]

Hi -- sorry for dropping the ball at the end of this PR. Thanks for merging it in! My hope is to implement a vector autoregression next. How would build the lag matrix without looping (theano.scan)? Some kind of algebraic operation? or some other theano reshape / broadcasting procedure?

[junpenglao]

@eph no problem ;-)
I have no given much thought about how to build the logp without theano.scan, but judging from the implementation from statsmodel it seems possible? Just build a large (sparse) covariance matrix?
Maybe we can leverage some of our GP modules?
We can further discuss it in #2491.