Hide keyboard shortcuts

Hot-keys on this page

r m x p   toggle line displays

j k   next/prev highlighted chunk

0   (zero) top of page

1   (one) first highlighted chunk

1""" 

2Innovations algorithm for MA(q) and SARIMA(p,d,q)x(P,D,Q,s) model parameters. 

3 

4Author: Chad Fulton 

5License: BSD-3 

6""" 

7import warnings 

8import numpy as np 

9 

10from scipy.optimize import minimize 

11from statsmodels.tools.tools import Bunch 

12from statsmodels.tsa.innovations import arma_innovations 

13from statsmodels.tsa.stattools import acovf, innovations_algo 

14from statsmodels.tsa.statespace.tools import diff 

15 

16from statsmodels.tsa.arima.specification import SARIMAXSpecification 

17from statsmodels.tsa.arima.params import SARIMAXParams 

18from statsmodels.tsa.arima.estimators.hannan_rissanen import hannan_rissanen 

19 

20 

21def innovations(endog, ma_order=0, demean=True): 

22 """ 

23 Estimate MA parameters using innovations algorithm. 

24 

25 Parameters 

26 ---------- 

27 endog : array_like or SARIMAXSpecification 

28 Input time series array, assumed to be stationary. 

29 ma_order : int, optional 

30 Maximum moving average order. Default is 0. 

31 demean : bool, optional 

32 Whether to estimate and remove the mean from the process prior to 

33 fitting the moving average coefficients. Default is True. 

34 

35 Returns 

36 ------- 

37 parameters : list of SARIMAXParams objects 

38 List elements correspond to estimates at different `ma_order`. For 

39 example, parameters[0] is an `SARIMAXParams` instance corresponding to 

40 `ma_order=0`. 

41 other_results : Bunch 

42 Includes one component, `spec`, containing the `SARIMAXSpecification` 

43 instance corresponding to the input arguments. 

44 

45 Notes 

46 ----- 

47 The primary reference is [1]_, section 5.1.3. 

48 

49 This procedure assumes that the series is stationary. 

50 

51 References 

52 ---------- 

53 .. [1] Brockwell, Peter J., and Richard A. Davis. 2016. 

54 Introduction to Time Series and Forecasting. Springer. 

55 """ 

56 max_spec = SARIMAXSpecification(endog, ma_order=ma_order) 

57 endog = max_spec.endog 

58 

59 if demean: 

60 endog = endog - endog.mean() 

61 

62 if not max_spec.is_ma_consecutive: 

63 raise ValueError('Innovations estimation unavailable for models with' 

64 ' seasonal or otherwise non-consecutive MA orders.') 

65 

66 sample_acovf = acovf(endog, fft=True) 

67 theta, v = innovations_algo(sample_acovf, nobs=max_spec.ma_order + 1) 

68 ma_params = [theta[i, :i] for i in range(1, max_spec.ma_order + 1)] 

69 sigma2 = v 

70 

71 out = [] 

72 for i in range(max_spec.ma_order + 1): 

73 spec = SARIMAXSpecification(ma_order=i) 

74 p = SARIMAXParams(spec=spec) 

75 if i == 0: 

76 p.params = sigma2[i] 

77 else: 

78 p.params = np.r_[ma_params[i - 1], sigma2[i]] 

79 out.append(p) 

80 

81 # Construct other results 

82 other_results = Bunch({ 

83 'spec': spec, 

84 }) 

85 

86 return out, other_results 

87 

88 

89def innovations_mle(endog, order=(0, 0, 0), seasonal_order=(0, 0, 0, 0), 

90 demean=True, enforce_invertibility=True, 

91 start_params=None, minimize_kwargs=None): 

92 """ 

93 Estimate SARIMA parameters by MLE using innovations algorithm. 

94 

95 Parameters 

96 ---------- 

97 endog : array_like 

98 Input time series array. 

99 order : tuple, optional 

100 The (p,d,q) order of the model for the number of AR parameters, 

101 differences, and MA parameters. Default is (0, 0, 0). 

102 seasonal_order : tuple, optional 

103 The (P,D,Q,s) order of the seasonal component of the model for the 

104 AR parameters, differences, MA parameters, and periodicity. Default 

105 is (0, 0, 0, 0). 

106 demean : bool, optional 

107 Whether to estimate and remove the mean from the process prior to 

108 fitting the SARIMA coefficients. Default is True. 

109 enforce_invertibility : bool, optional 

110 Whether or not to transform the MA parameters to enforce invertibility 

111 in the moving average component of the model. Default is True. 

112 start_params : array_like, optional 

113 Initial guess of the solution for the loglikelihood maximization. The 

114 AR polynomial must be stationary. If `enforce_invertibility=True` the 

115 MA poylnomial must be invertible. If not provided, default starting 

116 parameters are computed using the Hannan-Rissanen method. 

117 minimize_kwargs : dict, optional 

118 Arguments to pass to scipy.optimize.minimize. 

119 

120 Returns 

121 ------- 

122 parameters : SARIMAXParams object 

123 other_results : Bunch 

124 Includes four components: `spec`, containing the `SARIMAXSpecification` 

125 instance corresponding to the input arguments; `minimize_kwargs`, 

126 containing any keyword arguments passed to `minimize`; `start_params`, 

127 containing the untransformed starting parameters passed to `minimize`; 

128 and `minimize_results`, containing the output from `minimize`. 

129 

130 Notes 

131 ----- 

132 The primary reference is [1]_, section 5.2. 

133 

134 Note: we do not include `enforce_stationarity` as an argument, because this 

135 function requires stationarity. 

136 

137 TODO: support concentrating out the scale (should be easy: use sigma2=1 

138 and then compute sigma2=np.sum(u**2 / v) / len(u); would then need to 

139 redo llf computation in the Cython function). 

140 

141 TODO: add support for fixed parameters 

142 

143 TODO: add support for secondary optimization that does not enforce 

144 stationarity / invertibility, starting from first step's parameters 

145 

146 References 

147 ---------- 

148 .. [1] Brockwell, Peter J., and Richard A. Davis. 2016. 

149 Introduction to Time Series and Forecasting. Springer. 

150 """ 

151 spec = SARIMAXSpecification( 

152 endog, order=order, seasonal_order=seasonal_order, 

153 enforce_stationarity=True, enforce_invertibility=enforce_invertibility) 

154 endog = spec.endog 

155 if spec.is_integrated: 

156 warnings.warn('Provided `endog` series has been differenced to' 

157 ' eliminate integration prior to ARMA parameter' 

158 ' estimation.') 

159 endog = diff(endog, k_diff=spec.diff, 

160 k_seasonal_diff=spec.seasonal_diff, 

161 seasonal_periods=spec.seasonal_periods) 

162 if demean: 

163 endog = endog - endog.mean() 

164 

165 p = SARIMAXParams(spec=spec) 

166 

167 if start_params is None: 

168 sp = SARIMAXParams(spec=spec) 

169 

170 # Estimate starting parameters via Hannan-Rissanen 

171 hr, hr_results = hannan_rissanen(endog, ar_order=spec.ar_order, 

172 ma_order=spec.ma_order, demean=False) 

173 if spec.seasonal_periods == 0: 

174 # If no seasonal component, then `hr` gives starting parameters 

175 sp.params = hr.params 

176 else: 

177 # If we do have a seasonal component, estimate starting parameters 

178 # for the seasonal lags using the residuals from the previous step 

179 _ = SARIMAXSpecification( 

180 endog, seasonal_order=seasonal_order, 

181 enforce_stationarity=True, 

182 enforce_invertibility=enforce_invertibility) 

183 

184 ar_order = np.array(spec.seasonal_ar_lags) * spec.seasonal_periods 

185 ma_order = np.array(spec.seasonal_ma_lags) * spec.seasonal_periods 

186 seasonal_hr, seasonal_hr_results = hannan_rissanen( 

187 hr_results.resid, ar_order=ar_order, ma_order=ma_order, 

188 demean=False) 

189 

190 # Set the starting parameters 

191 sp.ar_params = hr.ar_params 

192 sp.ma_params = hr.ma_params 

193 sp.seasonal_ar_params = seasonal_hr.ar_params 

194 sp.seasonal_ma_params = seasonal_hr.ma_params 

195 sp.sigma2 = seasonal_hr.sigma2 

196 

197 # Then, require starting parameters to be stationary and invertible 

198 if not sp.is_stationary: 

199 sp.ar_params = [0] * sp.k_ar_params 

200 sp.seasonal_ar_params = [0] * sp.k_seasonal_ar_params 

201 

202 if not sp.is_invertible and spec.enforce_invertibility: 

203 sp.ma_params = [0] * sp.k_ma_params 

204 sp.seasonal_ma_params = [0] * sp.k_seasonal_ma_params 

205 

206 start_params = sp.params 

207 else: 

208 sp = SARIMAXParams(spec=spec) 

209 sp.params = start_params 

210 if not sp.is_stationary: 

211 raise ValueError('Given starting parameters imply a non-stationary' 

212 ' AR process. Innovations algorithm requires a' 

213 ' stationary process.') 

214 

215 if spec.enforce_invertibility and not sp.is_invertible: 

216 raise ValueError('Given starting parameters imply a non-invertible' 

217 ' MA process with `enforce_invertibility=True`.') 

218 

219 def obj(params): 

220 p.params = spec.constrain_params(params) 

221 

222 return -arma_innovations.arma_loglike( 

223 endog, ar_params=-p.reduced_ar_poly.coef[1:], 

224 ma_params=p.reduced_ma_poly.coef[1:], sigma2=p.sigma2) 

225 

226 # Untransform the starting parameters 

227 unconstrained_start_params = spec.unconstrain_params(start_params) 

228 

229 # Perform the minimization 

230 if minimize_kwargs is None: 

231 minimize_kwargs = {} 

232 if 'options' not in minimize_kwargs: 

233 minimize_kwargs['options'] = {} 

234 minimize_kwargs['options'].setdefault('maxiter', 100) 

235 minimize_results = minimize(obj, unconstrained_start_params, 

236 **minimize_kwargs) 

237 

238 # TODO: show warning if convergence failed. 

239 

240 # Reverse the transformation to get the optimal parameters 

241 p.params = spec.constrain_params(minimize_results.x) 

242 

243 # Construct other results 

244 other_results = Bunch({ 

245 'spec': spec, 

246 'minimize_results': minimize_results, 

247 'minimize_kwargs': minimize_kwargs, 

248 'start_params': start_params 

249 }) 

250 

251 return p, other_results