Yule-Walker estimation is simple, but it has some problems:
To obtain the Yule-Walker estimates for an AR model of order m, it is necessary to solve a system of linear equations with m unknowns. In addition, to select the order of the AR model by the minimum AIC method, we need to evaluate the AIC values of the models with orders up to M, the maximum order. Namely, we have to estimate the coefficients by solving systems of linear equations with one unknown, . . . , M unknowns.
So we have Durbin-Levinson Algorithm to deal with this problem. Hereinafter, the AR coefficients and the innovation variance of the AR model of order m are denoted as and
, respectively.The algorithm is defined as follows:
- Set
and
- For m = 1, …, M, repeat the following steps:
,
for i=1, …, m-1
,
Here we see a for the first, which means the
parameter in the order m. With this algorithm, we can easily estimate all parameters from order 1 to M in loop. The sample code is:
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 | long N = this.getSampleCount(); this.v[0] = getGamma(0); this.aic[0] = N * (Math.log(2 * Math.PI * this.v[0]) + 1) + 2; for (int m=1; m<=p; m++) { // estimate phi[m][m] double sum = 0; for (int j=1; j<=m-1;j++) { sum += this.phi[m-1][j] * getGamma(m-j); } this.phi[m][m] = (getGamma(m) - sum)/this.v[m-1]; // estimate phi[m][i], i=1,...,m-1 for (int i=1; i<=m-1; i++) { this.phi[m][i] = this.phi[m-1][i] - this.phi[m][m] * this.phi[m-1][m-i]; } // estimate v[m] this.v[m] = this.v[m-1] * (1 - this.phi[m][m] * this.phi[m][m]); this.aic[m] = N * (Math.log(2 * Math.PI * this.v[m]) + 1) + 2 * (m + 1); this.pacf[m] = this.phi[m][m]; } |
Because Durbin-Levinson Algorithm is tend to solve the same formula as Yule-Walker method, so it’s forecasting method is just the same as previous post.
The last implementation is Burg’s algorithm based on the maximum entropy method (MEM) to estimate. The process is:
- Set
, for n = 1, …, N. In addition, for the AR model of order 0, compute
, and
- For m = 1, …, M, repeat the following steps:
- Estimate the
by any of the formulae:
- Obtain the AR coefficients
, …,
by
- For n = m+1,…,N, obtain the forward prediction error as
- For n = m + 1, …, N, obtain the backward prediction error as
- Estimate the innovation variance of the AR model of order m by
- Obtain AIC by
- Estimate the
The following is sample code. Note here we use the third formula in step 2a.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 45 46 47 48 49 50 51 52 53 54 55 | double[] y = getSamples(); int N = y.length; double mean = this.getMean(); // let the mean to 0, double sumOfSquares = 0; double[] yy = new double[y.length]; for (int i=0;i<y.length;i++) { yy[i] = y[i] - mean; sumOfSquares += yy[i] * yy[i]; } double[][] vv = new double[this.p+1][y.length+1]; double[][] w = new double[this.p+1][y.length+1]; // initialize for (int i = 1; i<=yy.length; i++) { vv[0][i] = yy[i-1]; w[0][i] = yy[i-1]; } // estimate v[0] and aic[0] this.v[0] = sumOfSquares / N; //System.out.println(sumOfSquares); //System.out.println(this.getAcf().getStatistics().getSumsq()); //System.out.println(this.v[0]); this.aic[0] = N * (Math.log(2 * Math.PI * this.v[0] * this.v[0]) + 1) + 2; for (int m=1;m<=this.p; m++) { // estimate phi[m][m] double sum1 = 0; double sum2 = 0; double sum3 = 0; for (int n=m+1;n<=N;n++) { sum1 += vv[m-1][n] * w[m-1][n-m]; } for (int n=m+1;n<=N;n++) { sum2 += w[m-1][n-m] * w[m-1][n-m]; } for (int n=m+1;n<=N;n++) { sum3 += vv[m-1][n] * vv[m-1][n]; } this.phi[m][m] = 2 * (sum1 / (sum2 + sum3)); // estimate phi[m] .. phi[m][m-1] for (int j=1;j<=m-1;j++) { this.phi[m][j] = this.phi[m-1][j] - this.phi[m][m] * this.phi[m-1][m-j]; } // For n = m+1, ···,N, obtain the forward prediction for (int n=m+1; n<=N;n++) { vv[m][n] = vv[m-1][n] - this.phi[m][m] * w[m-1][n-m]; } // For n = m+1, ···,N, obtain the backward prediction error for (int n=m+1; n<=N;n++) { w[m][n-m] = w[m-1][n-m] - this.phi[m][m] * vv[m-1][n]; } // Estimate the innovation variance of the AR model of order m this.v[m] = this.v[m-1] * (1 - this.phi[m][m] * this.phi[m][m]); // Estimate aic this.aic[m] = N * (Math.log(2 * Math.PI * this.v[m]) + 1) + 2 * (m + 1); } |
Here is another sample implementation based on http://emptyloops.com/technotes/A%20tutorial%20on%20Burg’s%20method,%20algorithm%20and%20recursion.pdf.
1 2 3 4 5 6 7 8 9 10 11 12 13 14 15 16 17 18 19 20 21 22 23 24 25 26 27 28 29 30 31 32 33 34 35 36 37 38 39 40 41 42 43 44 | double[] x = getSamples(); int N = x.length -1; int m = this.p; // initialize Ak double[] Ak = new double[this.p+1]; Ak[0] = 1; // initialize f and b double[] f = ArrayUtils.clone(x); double[] b = ArrayUtils.clone(x); // initialize Dk double Dk = 0; for (int j=0;j<=N;j++) { Dk += 2 * f[j] * f[j]; } Dk -= f[0] * f[0] + b[N] * b[N]; // Burg Recursion for (int k=0; k< m; k++) { // compute mu double mu = 0; for (int n=0; n<=N-k-1;n++) { mu += f[n+k+1]*b[n]; } mu *= -2/Dk; // update Ak for (int n=0;n<=(k+1)/2;n++) { double t1 = Ak[n] + mu * Ak[k+1-n]; double t2 = Ak[k+1-n] + mu*Ak[n]; Ak[n] = t1; Ak[k+1-n] = t2; } // update f and b for (int n=0;n<=N-k-1;n++) { double t1 = f[n+k+1] + mu*b[n]; double t2 = b[n] + mu*f[n+k+1]; f[n+k+1] = t1; b[n] = t2; } // update Dk Dk = (1-mu*mu)*Dk - f[k+1]*f[k+1] - b[N-k-1]*b[N-k-1]; } //assign coefficients for (int i=1;i<=this.p;i++) { this.phi[this.p][i] = Ak[i] * -1; // remember to times -1 because this algorithm model is reversed } |
Note the sign of coefficients because the model is different.
Finally, we need to know the performance of them. In my simple simulation, the Burg’s algorithm has the best performance and Durbin-Levinson algorithm has the worst, but only a 1-2 ms difference in around 5,000 – 7,000 points. Yule-Walker method can only estimate an order at time, so if you need to test different orders to select the best one, it will cost another loop to estimate and might cost more time. As prediction, their 1 and 2 step forecasting is not significantly different.
Reference:
Genshiro Kitagawa, “Introduction to Time Series Modeling", http://www.amazon.com/Introduction-Modeling-Monographs-Statistics-Probability/dp/1584889217/ref=sr_1_1?ie=UTF8&qid=1361603954&sr=8-1&keywords=introduction+to+time+series+modeling
Related Free Lance Projects on GAF:
Start discussion »