2
Created on Sun Oct 10 14:57:50 2010
4
Author: josef-pktd, Skipper Seabold
7
TODO: check everywhere initialization of signal.lfilter
12
from scipy import signal, optimize
13
from scikits.statsmodels.model import LikelihoodModel, GenericLikelihoodModel
16
#copied from sandbox/regression/mle.py
17
#rename until merge of classes is complete
18
class Arma(GenericLikelihoodModel): #switch to generic mle
20
univariate Autoregressive Moving Average model
22
Note: This is not working yet, or does it
23
this can subclass TSMLEModel
26
def __init__(self, endog, exog=None):
27
#need to override p,q (nar,nma) correctly
28
super(Arma, self).__init__(endog, exog)
29
#set default arma(1,1)
38
def geterrors(self, params):
39
#copied from sandbox.tsa.arima.ARIMA
40
p, q = self.nar, self.nma
41
rhoy = np.concatenate(([1], params[:p]))
42
rhoe = np.concatenate(([1], params[p:p+q]))
43
errorsest = signal.lfilter(rhoy, rhoe, self.endog)
48
def loglike(self, params):
50
Loglikelihood for arma model
54
The ancillary parameter is assumed to be the last element of
58
# #copied from sandbox.tsa.arima.ARIMA
60
# rhoy = np.concatenate(([1], params[:p]))
61
# rhoe = np.concatenate(([1], params[p:-1]))
62
# errorsest = signal.lfilter(rhoy, rhoe, self.endog)
63
errorsest = self.geterrors(params)
64
sigma2 = np.maximum(params[-1]**2, 1e-6)
67
#this doesn't help for exploding paths
68
#errorsest[np.isnan(errorsest)] = 100
69
# llike = -0.5 * (np.sum(np.log(sigma2),axis)
70
# + np.sum((errorsest**2)/sigma2, axis)
71
# + nobs*np.log(2*np.pi))
72
llike = -0.5 * (nobs*np.log(sigma2)
73
+ np.sum((errorsest**2)/sigma2, axis)
74
+ nobs*np.log(2*np.pi))
77
#add for Jacobian calculation bsejac in GenericMLE, copied from loglike
78
def nloglikeobs(self, params):
80
Loglikelihood for arma model
84
The ancillary parameter is assumed to be the last element of
88
# #copied from sandbox.tsa.arima.ARIMA
90
# rhoy = np.concatenate(([1], params[:p]))
91
# rhoe = np.concatenate(([1], params[p:-1]))
92
# errorsest = signal.lfilter(rhoy, rhoe, self.endog)
93
errorsest = self.geterrors(params)
94
sigma2 = np.maximum(params[-1]**2, 1e-6)
97
#this doesn't help for exploding paths
98
#errorsest[np.isnan(errorsest)] = 100
99
# llike = -0.5 * (np.sum(np.log(sigma2),axis)
100
# + np.sum((errorsest**2)/sigma2, axis)
101
# + nobs*np.log(2*np.pi))
102
llike = 0.5 * (np.log(sigma2)
103
+ (errorsest**2)/sigma2
108
# def score(self, params):
110
# Score vector for Arma model
114
# jac = ndt.Jacobian(self.loglike, stepMax=1e-4)
115
# return jac(params)[-1]
119
# def hessian(self, params):
121
# Hessian of arma model. Currently uses numdifftools
124
# Hfun = ndt.Jacobian(self.score, stepMax=1e-4)
125
# return Hfun(params)[-1]
127
#copied from arima.ARIMA, needs splitting out of method specific code
128
def fit(self, order=(0,0,0), method="ls", rhoy0=None, rhoe0=None):
130
Estimate lag coefficients of an ARIMA process.
135
p,d,q where p is the number of AR lags, d is the number of
136
differences to induce stationarity, and q is the number of
138
method : str {"ls", "ssm"}
139
Method of estimation. LS is conditional least squares.
140
SSM is state-space model and the Kalman filter is used to
141
maximize the exact likelihood.
142
rhoy0, rhoe0 : array_like (optional)
143
starting values for estimation
147
rh, cov_x, infodict, mesg, ier : output of scipy.optimize.leastsq
149
estimate of lag parameters, concatenated [rhoy, rhoe]
151
unscaled (!) covariance matrix of coefficient estimates
153
if not hasattr(order, '__iter__'):
154
raise ValueError("order must be an iterable sequence. Got type \
155
%s instead" % type(order))
158
self.nar = p # needed for geterrors, needs cleanup
162
raise ValueError("Differencing not implemented yet")
163
# assume no constant, ie mu = 0
164
# unless overwritten then use w_bar for mu
165
Y = np.diff(endog, d, axis=0) #TODO: handle lags?
167
x = self.endog.squeeze() # remove the squeeze might be needed later
170
# rhoy = np.concatenate(([1], rho[:p]))
171
# rhoe = np.concatenate(([1], rho[p:]))
172
# etahatr = signal.lfilter(rhoy, rhoe, x)
173
# #print rho,np.sum(etahatr*etahatr)
176
#replace with start_params
178
rhoy0 = 0.5 * np.ones(p)
180
rhoe0 = 0.5 * np.ones(q)
182
method = method.lower()
185
#changes: use self.geterrors (nobs,):
186
# rh, cov_x, infodict, mesg, ier = \
187
# optimize.leastsq(errfn, np.r_[rhoy0, rhoe0],ftol=1e-10,full_output=True)
188
rh, cov_x, infodict, mesg, ier = \
189
optimize.leastsq(self.geterrors, np.r_[rhoy0, rhoe0],ftol=1e-10,full_output=True)
190
#TODO: need missing parameter estimates for LS, scale, residual-sdt
191
#TODO: integrate this into the MLE.fit framework?
192
elif method == "ssm":
194
else: #this is also conditional least squares
195
# fmin_bfgs is slow or doesn't work yet
196
errfnsum = lambda rho : np.sum(self.geterrors(rho)**2)
197
#xopt, {fopt, gopt, Hopt, func_calls, grad_calls
198
rh, fopt, gopt, cov_x, _,_, ier = \
199
optimize.fmin_bfgs(errfnsum, np.r_[rhoy0, rhoe0], maxiter=2, full_output=True)
200
infodict, mesg = None, None
202
self.ar_est = np.concatenate(([1], rh[:p]))
203
self.ma_est = np.concatenate(([1], rh[p:])) #rh[-q:])) doesnt work for q=0
204
self.error_estimate = self.geterrors(rh)
205
return rh, cov_x, infodict, mesg, ier
208
#renamed and needs check with other fit
209
def fit_mle(self, start_params=None, maxiter=5000, method='nm', tol=1e-08):
210
if start_params is None:
211
start_params = np.concatenate((0.05*np.ones(self.nar + self.nma), [1]))
212
mlefit = super(Arma, self).fit(start_params=start_params,
213
maxiter=maxiter, method=method, tol=tol)
216
#copied from arima.ARIMA
217
def predicted(self, ar=None, ma=None):
218
'''past predicted values of time series
219
just added, not checked yet
222
# #ar, ma not used, not useful as arguments for predicted pattern
223
# #need it for prediction for other time series, endog
228
return self.x + self.error_estimate
230
#copied from arima.ARIMA
231
def forecast(self, ar=None, ma=None, nperiod=10):
232
'''nperiod ahead forecast at the end of the data period
234
eta = np.r_[self.error_estimate, np.zeros(nperiod)]
239
return signal.lfilter(ma, ar, eta)
241
#copied from arima.ARIMA
242
#TODO: is this needed as a method at all?
243
#JP: not needed in this form, but can be replace with using the parameters
245
def generate_sample(cls, ar, ma, nsample, std=1):
246
eta = std * np.random.randn(nsample)
247
return signal.lfilter(ma, ar, eta)
added
removed
scikits/statsmodels/tsa/arma_mle.py
