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

1import numpy as np 

2 

3from statsmodels.tsa import arima_process 

4from statsmodels.tsa.statespace.tools import prefix_dtype_map 

5from statsmodels.tools.numdiff import _get_epsilon, approx_fprime_cs 

6from scipy.linalg.blas import find_best_blas_type 

7from . import _arma_innovations 

8 

9 

10def arma_innovations(endog, ar_params=None, ma_params=None, sigma2=1, 

11 normalize=False, prefix=None): 

12 """ 

13 Compute innovations using a given ARMA process. 

14 

15 Parameters 

16 ---------- 

17 endog : ndarray 

18 The observed time-series process, may be univariate or multivariate. 

19 ar_params : ndarray, optional 

20 Autoregressive parameters. 

21 ma_params : ndarray, optional 

22 Moving average parameters. 

23 sigma2 : ndarray, optional 

24 The ARMA innovation variance. Default is 1. 

25 normalize : bool, optional 

26 Whether or not to normalize the returned innovations. Default is False. 

27 prefix : str, optional 

28 The BLAS prefix associated with the datatype. Default is to find the 

29 best datatype based on given input. This argument is typically only 

30 used internally. 

31 

32 Returns 

33 ------- 

34 innovations : ndarray 

35 Innovations (one-step-ahead prediction errors) for the given `endog` 

36 series with predictions based on the given ARMA process. If 

37 `normalize=True`, then the returned innovations have been "whitened" by 

38 dividing through by the square root of the mean square error. 

39 innovations_mse : ndarray 

40 Mean square error for the innovations. 

41 """ 

42 # Parameters 

43 endog = np.array(endog) 

44 squeezed = endog.ndim == 1 

45 if squeezed: 

46 endog = endog[:, None] 

47 

48 ar_params = np.atleast_1d([] if ar_params is None else ar_params) 

49 ma_params = np.atleast_1d([] if ma_params is None else ma_params) 

50 

51 nobs, k_endog = endog.shape 

52 ar = np.r_[1, -ar_params] 

53 ma = np.r_[1, ma_params] 

54 

55 # Get BLAS prefix 

56 if prefix is None: 

57 prefix, dtype, _ = find_best_blas_type( 

58 [endog, ar_params, ma_params, np.array(sigma2)]) 

59 dtype = prefix_dtype_map[prefix] 

60 

61 # Make arrays contiguous for BLAS calls 

62 endog = np.asfortranarray(endog, dtype=dtype) 

63 ar_params = np.asfortranarray(ar_params, dtype=dtype) 

64 ma_params = np.asfortranarray(ma_params, dtype=dtype) 

65 sigma2 = dtype(sigma2).item() 

66 

67 # Get the appropriate functions 

68 arma_transformed_acovf_fast = getattr( 

69 _arma_innovations, prefix + 'arma_transformed_acovf_fast') 

70 arma_innovations_algo_fast = getattr( 

71 _arma_innovations, prefix + 'arma_innovations_algo_fast') 

72 arma_innovations_filter = getattr( 

73 _arma_innovations, prefix + 'arma_innovations_filter') 

74 

75 # Run the innovations algorithm for ARMA coefficients 

76 arma_acovf = arima_process.arma_acovf(ar, ma, 

77 sigma2=sigma2, nobs=nobs) / sigma2 

78 acovf, acovf2 = arma_transformed_acovf_fast(ar, ma, arma_acovf) 

79 theta, v = arma_innovations_algo_fast(nobs, ar_params, ma_params, 

80 acovf, acovf2) 

81 v = np.array(v) 

82 if normalize: 

83 v05 = v**0.5 

84 

85 # Run the innovations filter across each series 

86 u = [] 

87 for i in range(k_endog): 

88 u_i = np.array(arma_innovations_filter(endog[:, i], ar_params, 

89 ma_params, theta)) 

90 u.append(u_i / v05 if normalize else u_i) 

91 u = np.vstack(u).T 

92 

93 # Post-processing 

94 if squeezed: 

95 u = u.squeeze() 

96 

97 return u, v 

98 

99 

100def arma_loglike(endog, ar_params=None, ma_params=None, sigma2=1, prefix=None): 

101 """ 

102 Compute the log-likelihood of the given data assuming an ARMA process. 

103 

104 Parameters 

105 ---------- 

106 endog : ndarray 

107 The observed time-series process. 

108 ar_params : ndarray, optional 

109 Autoregressive parameters. 

110 ma_params : ndarray, optional 

111 Moving average parameters. 

112 sigma2 : ndarray, optional 

113 The ARMA innovation variance. Default is 1. 

114 prefix : str, optional 

115 The BLAS prefix associated with the datatype. Default is to find the 

116 best datatype based on given input. This argument is typically only 

117 used internally. 

118 

119 Returns 

120 ------- 

121 float 

122 The joint loglikelihood. 

123 """ 

124 llf_obs = arma_loglikeobs(endog, ar_params=ar_params, ma_params=ma_params, 

125 sigma2=sigma2, prefix=prefix) 

126 return np.sum(llf_obs) 

127 

128 

129def arma_loglikeobs(endog, ar_params=None, ma_params=None, sigma2=1, 

130 prefix=None): 

131 """ 

132 Compute the log-likelihood for each observation assuming an ARMA process. 

133 

134 Parameters 

135 ---------- 

136 endog : ndarray 

137 The observed time-series process. 

138 ar_params : ndarray, optional 

139 Autoregressive parameters. 

140 ma_params : ndarray, optional 

141 Moving average parameters. 

142 sigma2 : ndarray, optional 

143 The ARMA innovation variance. Default is 1. 

144 prefix : str, optional 

145 The BLAS prefix associated with the datatype. Default is to find the 

146 best datatype based on given input. This argument is typically only 

147 used internally. 

148 

149 Returns 

150 ------- 

151 ndarray 

152 Array of loglikelihood values for each observation. 

153 """ 

154 endog = np.array(endog) 

155 ar_params = np.atleast_1d([] if ar_params is None else ar_params) 

156 ma_params = np.atleast_1d([] if ma_params is None else ma_params) 

157 

158 if prefix is None: 

159 prefix, dtype, _ = find_best_blas_type( 

160 [endog, ar_params, ma_params, np.array(sigma2)]) 

161 dtype = prefix_dtype_map[prefix] 

162 

163 endog = np.ascontiguousarray(endog, dtype=dtype) 

164 ar_params = np.asfortranarray(ar_params, dtype=dtype) 

165 ma_params = np.asfortranarray(ma_params, dtype=dtype) 

166 sigma2 = dtype(sigma2).item() 

167 

168 func = getattr(_arma_innovations, prefix + 'arma_loglikeobs_fast') 

169 return func(endog, ar_params, ma_params, sigma2) 

170 

171 

172def arma_score(endog, ar_params=None, ma_params=None, sigma2=1, 

173 prefix=None): 

174 """ 

175 Compute the score (gradient of the log-likelihood function). 

176 

177 Parameters 

178 ---------- 

179 endog : ndarray 

180 The observed time-series process. 

181 ar_params : ndarray, optional 

182 Autoregressive coefficients, not including the zero lag. 

183 ma_params : ndarray, optional 

184 Moving average coefficients, not including the zero lag, where the sign 

185 convention assumes the coefficients are part of the lag polynomial on 

186 the right-hand-side of the ARMA definition (i.e. they have the same 

187 sign from the usual econometrics convention in which the coefficients 

188 are on the right-hand-side of the ARMA definition). 

189 sigma2 : ndarray, optional 

190 The ARMA innovation variance. Default is 1. 

191 prefix : str, optional 

192 The BLAS prefix associated with the datatype. Default is to find the 

193 best datatype based on given input. This argument is typically only 

194 used internally. 

195 

196 Returns 

197 --------- 

198 ndarray 

199 Score, evaluated at the given parameters. 

200 

201 Notes 

202 ----- 

203 This is a numerical approximation, calculated using first-order complex 

204 step differentiation on the `arma_loglike` method. 

205 """ 

206 ar_params = [] if ar_params is None else ar_params 

207 ma_params = [] if ma_params is None else ma_params 

208 

209 p = len(ar_params) 

210 q = len(ma_params) 

211 

212 def func(params): 

213 return arma_loglike(endog, params[:p], params[p:p + q], params[p + q:]) 

214 

215 params0 = np.r_[ar_params, ma_params, sigma2] 

216 epsilon = _get_epsilon(params0, 2., None, len(params0)) 

217 return approx_fprime_cs(params0, func, epsilon) 

218 

219 

220def arma_scoreobs(endog, ar_params=None, ma_params=None, sigma2=1, 

221 prefix=None): 

222 """ 

223 Compute the score (gradient) per observation. 

224 

225 Parameters 

226 ---------- 

227 endog : ndarray 

228 The observed time-series process. 

229 ar_params : ndarray, optional 

230 Autoregressive coefficients, not including the zero lag. 

231 ma_params : ndarray, optional 

232 Moving average coefficients, not including the zero lag, where the sign 

233 convention assumes the coefficients are part of the lag polynomial on 

234 the right-hand-side of the ARMA definition (i.e. they have the same 

235 sign from the usual econometrics convention in which the coefficients 

236 are on the right-hand-side of the ARMA definition). 

237 sigma2 : ndarray, optional 

238 The ARMA innovation variance. Default is 1. 

239 prefix : str, optional 

240 The BLAS prefix associated with the datatype. Default is to find the 

241 best datatype based on given input. This argument is typically only 

242 used internally. 

243 

244 Returns 

245 --------- 

246 ndarray 

247 Score per observation, evaluated at the given parameters. 

248 

249 Notes 

250 ----- 

251 This is a numerical approximation, calculated using first-order complex 

252 step differentiation on the `arma_loglike` method. 

253 """ 

254 ar_params = [] if ar_params is None else ar_params 

255 ma_params = [] if ma_params is None else ma_params 

256 

257 p = len(ar_params) 

258 q = len(ma_params) 

259 

260 def func(params): 

261 return arma_loglikeobs(endog, params[:p], params[p:p + q], 

262 params[p + q:]) 

263 

264 params0 = np.r_[ar_params, ma_params, sigma2] 

265 epsilon = _get_epsilon(params0, 2., None, len(params0)) 

266 return approx_fprime_cs(params0, func, epsilon)