pythonarma模型t检验_python时间序列分析（ARIMA模型）

1 #Note: The information criteria add 1 to the number of parameters

2 #whenever the model has an AR or MA term since, in principle,

3 #the variance could be treated as a free parameter and restricted

4 #This code does not allow this, but it adds consistency with other

5 #packages such as gretl and X12-ARIMA

6

7 from __future__ importabsolute_import8 from statsmodels.compat.python importstring_types, range9 #for 2to3 with extensions

10

11 from datetime importdatetime12

13 importnumpy as np14 from scipy importoptimize15 from scipy.stats importt, norm16 from scipy.signal importlfilter17 from numpy importdot, log, zeros, pi18 from numpy.linalg importinv19

20 from statsmodels.tools.decorators import(cache_readonly,21 resettable_cache)22 importstatsmodels.tsa.base.tsa_model as tsbase23 importstatsmodels.base.wrapper as wrap24 from statsmodels.regression.linear_model importyule_walker, GLS25 from statsmodels.tsa.tsatools import(lagmat, add_trend,26 _ar_transparams, _ar_invtransparams,27 _ma_transparams, _ma_invtransparams,28 unintegrate, unintegrate_levels)29 from statsmodels.tsa.vector_ar importutil30 from statsmodels.tsa.ar_model importAR31 from statsmodels.tsa.arima_process importarma2ma32 from statsmodels.tools.numdiff importapprox_hess_cs, approx_fprime_cs33 from statsmodels.tsa.base.datetools import_index_date34 from statsmodels.tsa.kalmanf importKalmanFilter35

36 _armax_notes = """

37

38 Notes39 -----40 If exogenous variables are given, then the model that is fit is41

42 .. math::43

44 \\phi(L)(y_t - X_t\\beta) = \\theta(L)\epsilon_t45

46 where :math:`\\phi` and :math:`\\theta` are polynomials in the lag47 operator, :math:`L`. This is the regression model with ARMA errors,48 or ARMAX model. This specification is used, whether or not the model49 is fit using conditional sum of square or maximum-likelihood, using50 the `method` argument in51 :meth:`statsmodels.tsa.arima_model.%(Model)s.fit`. Therefore, for52 now, `css` and `mle` refer to estimation methods only. This may53 change for the case of the `css` model in future versions.54 """

55

56 _arma_params = """\57 endog : array-like58 The endogenous variable.59 order : iterable60 The (p,q) order of the model for the number of AR parameters,61 differences, and MA parameters to use.62 exog : array-like, optional63 An optional arry of exogenous variables. This should *not* include a64 constant or trend. You can specify this in the `fit` method."""

65

66 _arma_model = "Autoregressive Moving Average ARMA(p,q) Model"

67

68 _arima_model = "Autoregressive Integrated Moving Average ARIMA(p,d,q) Model"

69

70 _arima_params = """\71 endog : array-like72 The endogenous variable.73 order : iterable74 The (p,d,q) order of the model for the number of AR parameters,75 differences, and MA parameters to use.76 exog : array-like, optional77 An optional arry of exogenous variables. This should *not* include a78 constant or trend. You can specify this in the `fit` method."""

79

80 _predict_notes = """

81 Notes82 -----83 Use the results predict method instead.84 """

85

86 _results_notes = """

87 Notes88 -----89 It is recommended to use dates with the time-series models, as the90 below will probably make clear. However, if ARIMA is used without91 dates and/or `start` and `end` are given as indices, then these92 indices are in terms of the *original*, undifferenced series. Ie.,93 given some undifferenced observations::94

95 1970Q1, 196 1970Q2, 1.597 1970Q3, 1.2598 1970Q4, 2.2599 1971Q1, 1.2100 1971Q2, 4.1101

102 1970Q1 is observation 0 in the original series. However, if we fit an103 ARIMA(p,1,q) model then we lose this first observation through104 differencing. Therefore, the first observation we can forecast (if105 using exact MLE) is index 1. In the differenced series this is index106 0, but we refer to it as 1 from the original series.107 """

108

109 _predict = """

110 %(Model)s model in-sample and out-of-sample prediction111

112 Parameters113 ----------114 %(params)s115 start : int, str, or datetime116 Zero-indexed observation number at which to start forecasting, ie.,117 the first forecast is start. Can also be a date string to118 parse or a datetime type.119 end : int, str, or datetime120 Zero-indexed observation number at which to end forecasting, ie.,121 the first forecast is start. Can also be a date string to122 parse or a datetime type. However, if the dates index does not123 have a fixed frequency, end must be an integer index if you124 want out of sample prediction.125 exog : array-like, optional126 If the model is an ARMAX and out-of-sample forecasting is127 requested, exog must be given. Note that you'll need to pass128 `k_ar` additional lags for any exogenous variables. E.g., if you129 fit an ARMAX(2, q) model and want to predict 5 steps, you need 7130 observations to do this.131 dynamic : bool, optional132 The `dynamic` keyword affects in-sample prediction. If dynamic133 is False, then the in-sample lagged values are used for134 prediction. If `dynamic` is True, then in-sample forecasts are135 used in place of lagged dependent variables. The first forecasted136 value is `start`.137 %(extra_params)s138

139 Returns140 -------141 %(returns)s142 %(extra_section)s143 """

144

145 _predict_returns = """predict : array146 The predicted values.147

148 """

149

150 _arma_predict = _predict % {"Model" : "ARMA",151 "params" : """

152 params : array-like153 The fitted parameters of the model.""",154 "extra_params" : "",155 "returns": _predict_returns,156 "extra_section": _predict_notes}157

158 _arma_results_predict = _predict % {"Model" : "ARMA", "params" : "",159 "extra_params" : "",160 "returns": _predict_returns,161 "extra_section": _results_notes}162

163 _arima_predict = _predict % {"Model" : "ARIMA",164 "params" : """params : array-like165 The fitted parameters of the model.""",166 "extra_params" : """typ : str {'linear', 'levels'}167

168 - 'linear' : Linear prediction in terms of the differenced169 endogenous variables.170 - 'levels' : Predict the levels of the original endogenous171 variables.\n""", "returns": _predict_returns,172 "extra_section": _predict_notes}173

174 _arima_results_predict = _predict % {"Model" : "ARIMA",175 "params" : "",176 "extra_params":177 """typ : str {'linear', 'levels'}178

179 - 'linear' : Linear prediction in terms of the differenced180 endogenous variables.181 - 'levels' : Predict the levels of the original endogenous182 variables.\n""",183 "returns": _predict_returns,184 "extra_section": _results_notes}185

186 _arima_plot_predict_example = """Examples187 --------188 >>> import statsmodels.api as sm189 >>> import matplotlib.pyplot as plt190 >>> import pandas as pd191 >>>192 >>> dta = sm.datasets.sunspots.load_pandas().data[['SUNACTIVITY']]193 >>> dta.index = pd.DatetimeIndex(start='1700', end='2009', freq='A')194 >>> res = sm.tsa.ARMA(dta, (3, 0)).fit()195 >>> fig, ax = plt.subplots()196 >>> ax = dta.ix['1950':].plot(ax=ax)197 >>> fig = res.plot_predict('1990', '2012', dynamic=True, ax=ax,198 ... plot_insample=False)199 >>> plt.show()200

201 .. plot:: plots/arma_predict_plot.py202 """

203

204 _plot_predict = ("""

205 Plot forecasts206 """ + '\n'.join(_predict.split('\n')[2:])) %{207 "params" : "",208 "extra_params" : """alpha : float, optional209 The confidence intervals for the forecasts are (1 - alpha)%210 plot_insample : bool, optional211 Whether to plot the in-sample series. Default is True.212 ax : matplotlib.Axes, optional213 Existing axes to plot with.""",214 "returns" : """fig : matplotlib.Figure215 The plotted Figure instance""",216 "extra_section" : ('\n' + _arima_plot_predict_example +

217 '\n' +_results_notes)218 }219

220 _arima_plot_predict = ("""

221 Plot forecasts222 """ + '\n'.join(_predict.split('\n')[2:])) %{223 "params" : "",224 "extra_params" : """alpha : float, optional225 The confidence intervals for the forecasts are (1 - alpha)%226 plot_insample : bool, optional227 Whether to plot the in-sample series. Default is True.228 ax : matplotlib.Axes, optional229 Existing axes to plot with.""",230 "returns" : """fig : matplotlib.Figure231 The plotted Figure instance""",232 "extra_section" : ('\n' + _arima_plot_predict_example +

233 '\n' +

234 '\n'.join(_results_notes.split('\n')[:3]) +

235 ("""

236 This is hard-coded to only allow plotting of the forecasts in levels.237 """) +

238 '\n'.join(_results_notes.split('\n')[3:]))239 }240

241

242 defcumsum_n(x, n):243 ifn:244 n -= 1

245 x =np.cumsum(x)246 returncumsum_n(x, n)247 else:248 returnx249

250

251 def_check_arima_start(start, k_ar, k_diff, method, dynamic):252 if start <0:253 raise ValueError("The start index %d of the original series"

254 "has been differenced away" %start)255 elif (dynamic or 'mle' not in method) and start <256 raise valueerror must be>= k_ar for conditional MLE"256>

257 "or dynamic forecast. Got %d" %start)258

259

260 def_get_predict_out_of_sample(endog, p, q, k_trend, k_exog, start, errors,261 trendparam, exparams, arparams, maparams, steps,262 method, exog=None):263 """

264 Returns endog, resid, mu of appropriate length for out of sample265 prediction.266 """

267 ifq:268 resid =np.zeros(q)269 if start and 'mle' in method or (start == p and not start ==0):270 resid[:q] = errors[start-q:start]271 elifstart:272 resid[:q] = errors[start-q-p:start-p]273 else:274 resid[:q] = errors[-q:]275 else:276 resid =None277

278 y =endog279 if k_trend == 1:280 #use expectation not constant

281 if k_exog >0:282 #TODO: technically should only hold for MLE not

283 #conditional model. See #274.

284 #ensure 2-d for conformability

285 if np.ndim(exog) == 1 and k_exog == 1:286 #have a 1d series of observations -> 2d

287 exog =exog[:, None]288 elif np.ndim(exog) == 1:289 #should have a 1d row of exog -> 2d

290 if len(exog) !=k_exog:291 raise ValueError("1d exog given and len(exog) != k_exog")292 exog =exog[None, :]293 X = lagmat(np.dot(exog, exparams), p, original='in', trim='both')294 mu = trendparam * (1 -arparams.sum())295 #arparams were reversed in unpack for ease later

296 mu = mu + (np.r_[1, -arparams[::-1]] * X).sum(1)[:, None]297 else:298 mu = trendparam * (1 -arparams.sum())299 mu = np.array([mu]*steps)300 elif k_exog >0:301 X =np.dot(exog, exparams)302 #NOTE: you shouldn't have to give in-sample exog!

303 X = lagmat(X, p, original='in', trim='both')304 mu = (np.r_[1, -arparams[::-1]] * X).sum(1)[:, None]305 else:306 mu =np.zeros(steps)307

308 endog = np.zeros(p + steps - 1)309

310 if p andstart:311 endog[:p] = y[start-p:start]312 elifp:313 endog[:p] = y[-p:]314

315 returnendog, resid, mu316

317

318 def_arma_predict_out_of_sample(params, steps, errors, p, q, k_trend, k_exog,319 endog, exog=None, start=0, method='mle'):320 (trendparam, exparams,321 arparams, maparams) =_unpack_params(params, (p, q), k_trend,322 k_exog, reverse=True)323 #print 'params:',params

324 #print 'arparams:',arparams,'maparams:',maparams

325 endog, resid, mu =_get_predict_out_of_sample(endog, p, q, k_trend, k_exog,326 start, errors, trendparam,327 exparams, arparams,328 maparams, steps, method,329 exog)330 #print 'mu[-1]:',mu[-1], 'mu[0]:',mu[0]

331 forecast =np.zeros(steps)332 if steps == 1:333 ifq:334 return mu[0] + np.dot(arparams, endog[:p]) +np.dot(maparams,335 resid[:q]), mu[0]336 else:337 return mu[0] +np.dot(arparams, endog[:p]), mu[0]338

339 ifq:340 i = 0 #if q == 1

341 else:342 i = -1

343

344 for i in range(min(q, steps - 1)):345 fcast = (mu[i] + np.dot(arparams, endog[i:i + p]) +

346 np.dot(maparams[:q - i], resid[i:i +q]))347 forecast[i] =fcast348 endog[i+p] =fcast349

350 for i in range(i + 1, steps - 1):351 fcast = mu[i] + np.dot(arparams, endog[i:i+p])352 forecast[i] =fcast353 endog[i+p] =fcast354

355 #need to do one more without updating endog

356 forecast[-1] = mu[-1] + np.dot(arparams, endog[steps - 1:])357 return forecast, mu[-1] #Modified by me, the former is return forecast

358

359

360 def_arma_predict_in_sample(start, end, endog, resid, k_ar, method):361 """

362 Pre- and in-sample fitting for ARMA.363 """

364 if 'mle' inmethod:365 fittedvalues = endog - resid #get them all then trim

366 else:367 fittedvalues = endog[k_ar:] -resid368

369 fv_start =start370 if 'mle' not inmethod:371 fv_start -= k_ar #start is in terms of endog index

372 fv_end = min(len(fittedvalues), end + 1)373 returnfittedvalues[fv_start:fv_end]374

375

376 def_validate(start, k_ar, k_diff, dates, method):377