1 static int tw_acf (const double *phi, int p,
2 const double *theta, int q,
3 double *acf, double *cvli, int ma,
4 double *alpha, int mxpq);
5
6 #define max0(i,j) ((i)>(j) ? (i) : (j))
7 #define min0(i,j) ((i)<(j) ? (i) : (j))
8
9 #define DSHRINK 0.0625
10 #define DEXPAND 16.0
11 #define DCDELTA 4.0
12
13 /*
14 Algorithm AS 197 Applied Statistics (1984) Vol. 33, No. 1
15
16 Computes a quantity monotonically related to the likelihood
17 function of an autoregressive-moving average process,
18 expressed as fact * sumsq.
19
20 Sign of the MA coefficients switched to agree with the
21 convention followed by gretl; notation revised somewhat
22 to be closer to Melard's documentation.
23 */
24
flikam(const double * phi,int p,const double * theta,int q,const double * w,double * e,int n,double * sumsq,double * fact,double * vw,double * vl,int rp1,double * vk,int r,double toler)25 int flikam (const double *phi, int p,
26 const double *theta, int q,
27 const double *w, double *e, int n,
28 double *sumsq, double *fact,
29 double *vw, double *vl, int rp1,
30 double *vk, int r, double toler)
31 {
32 double eps1 = 1.0e-10;
33 double detman = 1.0;
34 double detcar = 0.0;
35 double h2, a, aor, wi;
36 double vw0, vl0, alfa, flj;
37 int mxpq = max0(p, q);
38 int mnpq = min0(p, q);
39 int mxpqp1 = mxpq + 1;
40 int do_quick = 0;
41 int last, loop, jfrom, nexti;
42 int i, j, k, ret = 0;
43
44 *sumsq = *fact = 0;
45
46 /* calculation of the autocovariance function of a process with
47 unit innovation variance (vw) and the covariance between the
48 variable and the lagged innovations (vl).
49 */
50 ret = tw_acf(phi, p, theta, q, vw, vl, mxpqp1, vk, mxpq);
51 if (ret > 0) {
52 return ret;
53 }
54
55 /* computation of vk */
56 vk[0] = vw[0];
57 for (k=1; k<r; k++) {
58 vk[k] = 0;
59 for (j=k; j<p; j++) {
60 vk[k] += phi[j] * vw[j+1-k];
61 }
62 for (j=k; j<=q; j++) {
63 vk[k] += theta[j-1] * vl[j-k];
64 }
65 }
66
67 /* computation of the initial vectors (vl, vk) */
68 h2 = vk[0];
69 if (h2 <= eps1) {
70 return 8; /* fault code */
71 }
72 vl[r-1] = 0;
73 for (j=0; j<r; j++) {
74 vw[j] = 0;
75 if (j < r-1) {
76 vl[j] = vk[j+1];
77 }
78 if (j < p) {
79 vl[j] += phi[j] * h2;
80 }
81 vk[j] = vl[j];
82 }
83
84 /* exit if no observations */
85 if (n <= 0) {
86 return 9; /* fault code */
87 }
88
89 /* initialization for loop */
90 last = p - q;
91 loop = p;
92 jfrom = p;
93 vw[p] = 0;
94 vl[mxpq] = 0;
95
96 /* loop on time */
97 for (i=0; i<n; i++) {
98 /* test for skipping updating */
99 if (i == last) {
100 loop = mnpq;
101 jfrom = loop;
102 /* test for switching */
103 if (q == 0) {
104 /* pure AR */
105 do_quick = 1;
106 break;
107 }
108 }
109 if (i > mxpq && fabs(h2 - 1.0) < toler) {
110 do_quick = 1;
111 break;
112 }
113 /* update scalars */
114 detman *= h2;
115 while (fabs(detman) >= 1.0) {
116 detman *= DSHRINK;
117 detcar += DCDELTA;
118 }
119 while (fabs(detman) < DSHRINK) {
120 detman *= DEXPAND;
121 detcar -= DCDELTA;
122 }
123 vw0 = vw[0];
124 wi = w[i];
125 a = wi - vw0;
126 e[i] = a / sqrt(h2);
127 aor = a / h2;
128 *sumsq += a * aor;
129 vl0 = vl[0];
130 alfa = vl0 / h2;
131 h2 -= alfa * vl0;
132 if (h2 <= eps1) {
133 return 8; /* fault code */
134 }
135 /* update vectors */
136 for (j=0; j<loop; j++) {
137 flj = vl[j+1] + phi[j] * vl0;
138 vw[j] = vw[j+1] + phi[j] * vw0 + aor * vk[j];
139 vl[j] = flj - alfa * vk[j];
140 vk[j] -= alfa * flj;
141 }
142 for (j=jfrom; j<q; j++) {
143 vw[j] = vw[j+1] + aor * vk[j];
144 vl[j] = vl[j+1] - alfa * vk[j];
145 vk[j] -= alfa * vl[j+1];
146 }
147 for (j=jfrom; j<p; j++) {
148 vw[j] = vw[j+1] + phi[j] * wi;
149 }
150 }
151
152 if (do_quick) {
153 /* quick recursions */
154 nexti = i;
155 ret = -nexti;
156 for (i=nexti; i<n; i++) {
157 e[i] = w[i];
158 for (j=0; j<p; j++) {
159 e[i] -= phi[j] * w[i-j-1];
160 }
161 for (j=0; j<q; j++) {
162 e[i] -= theta[j] * e[i-j-1];
163 }
164 *sumsq += e[i] * e[i];
165 }
166 }
167
168 *fact = pow(detman, 1.0 / n) * pow(2.0, detcar / n);
169
170 return ret;
171 }
172
173 /*
174 Algorithm AS 197.1 Applied Statistics (1984) VOL. 33, NO. 1
175
176 Implementation of the algorithm of G. Tunnicliffe-Wilson
177 (J. Statist. Comput. Simul. 8, 1979, 301-309) for computation of the
178 autocovariance function of an ARMA process of order (p, q) and unit
179 innovation variance. The autoregressive and moving average
180 coefficients are stored in vectors phi and theta, using Box and
181 Jenkins notation. On output, vector cvli contains the covariances
182 between the variable and the (K-1)-lagged innovation for K = 1, ...,
183 q + 1.
184
185 Dimensions: phi(p), theta(q), acf(ma), cvli(mxpqp1), alpha(mxpq)
186 */
187
tw_acf(const double * phi,int p,const double * theta,int q,double * acf,double * cvli,int ma,double * alpha,int mxpq)188 static int tw_acf (const double *phi, int p,
189 const double *theta, int q,
190 double *acf, double *cvli, int ma,
191 double *alpha, int mxpq)
192 {
193 double ak, div, eps2 = 1.0e-10;
194 int i, j, k, kc, j1;
195 int kp1, mikp, miim1p;
196
197 /* Initialization, and return if p = q = 0 */
198 acf[0] = cvli[0] = 1.0;
199 if (ma == 1) {
200 return 0;
201 }
202 for (i=1; i<ma; i++) {
203 acf[i] = cvli[i] = 0;
204 }
205 for (k=0; k<mxpq; k++) {
206 alpha[k] = 0;
207 }
208
209 /* computation of the ACF of the moving average part,
210 stored in @acf
211 */
212 for (k=0; k<q; k++) {
213 cvli[k+1] = acf[k+1] = theta[k];
214 kc = q - k - 1;
215 for (j=0; j<kc; j++) {
216 acf[k+1] += theta[j] * theta[j+k+1];
217 }
218 acf[0] += theta[k] * theta[k];
219 }
220
221 if (p == 0) {
222 /* no AR part */
223 return 0;
224 }
225
226 /* initialization of cvli */
227 for (k=0; k<p; k++) {
228 alpha[k] = cvli[k] = phi[k];
229 }
230
231 /* Tunnicliffe-Wilson's alpha and delta (delta stored in ACF
232 which is gradually overwritten)
233 */
234 for (k=0; k<mxpq; k++) {
235 kc = mxpq - k - 1;
236 if (kc < p) {
237 ak = alpha[kc];
238 div = 1.0 - ak * ak;
239 if (div < eps2) {
240 return 5;
241 }
242 if (kc == 0) {
243 break;
244 }
245 for (j=0; j<kc; j++) {
246 alpha[j] = (cvli[j] + ak * cvli[kc-j-1]) / div;
247 }
248 for (j=0; j<kc; j++) {
249 cvli[j] = alpha[j];
250 }
251 }
252 if (kc < q) {
253 j1 = max0(kc - p, 1);
254 for (j=j1; j<=kc; j++) {
255 acf[j] += acf[kc+1] * alpha[kc-j];
256 }
257 }
258 }
259
260 /* computation of Tunnicliffe-Wilson's nu (stored in cvli,
261 copied into acf)
262 */
263 acf[0] *= 0.5;
264 for (k=0; k<p; k++) {
265 ak = alpha[k];
266 kp1 = k + 1;
267 div = 1.0 - ak * ak;
268 for (j=0; j<=kp1; j++) {
269 cvli[j] = (acf[j] + ak * acf[k+1-j]) / div;
270 }
271 for (j=0; j<=kp1; j++) {
272 acf[j] = cvli[j];
273 }
274 }
275
276 /* Computation of acf (which is gradually overwritten) */
277 for (i=0; i<ma; i++) {
278 miim1p = min0(i, p);
279 for (j=0; j<miim1p; j++) {
280 acf[i] += phi[j] * acf[i-j-1];
281 }
282 }
283 acf[0] *= 2.0;
284
285 /* computation of cvli */
286 cvli[0] = 1.0;
287 for (k=0; k<q; k++) {
288 cvli[k+1] = theta[k];
289 mikp = min0(k+1, p);
290 for (j=0; j<mikp; j++) {
291 cvli[k+1] += phi[j] * cvli[k-j];
292 }
293 }
294
295 return 0;
296 }
297