package com.imsl.test.example.stat; import java.text.*; import com.imsl.stat.*; import com.imsl.math.*; /** *

* Fits an \(\text{ARMA}(2,1)\) to the Wolfer * sunspot data using the method of maximum likelihood.

*

* Consider the Wolfer Sunspot Data (Anderson 1971, p. 660) consisting of the * number of sunspots observed each year from 1749 through 1924. The data set * for this example consists of the number of sunspots observed from 1770 * through 1869. The \(ARMA(2, 1)\) model can be expressed as

*

* $$z_t = \theta _0 + \phi _1 z_{t - 1} + \phi _2 z_{t - 2} - \theta _1 A_{t - * 1} + A_t$$

*

* where the errors \(A_t\) are independently normally distributed with mean * zero and variance \(\sigma _A^2\).

*

* The maximum likelihood estimates *

* $${\rm{\hat \theta }}_{\rm{0}} {\rm{,\hat \phi }}_{\rm{1}} {\rm{,\hat \phi * }}_{\rm{2}} {\rm{, and \, \hat \theta }}_{\rm{1}}$$

* from ARMAMaxLikelihood are compared to the same estimates using * the method of moments and least squares from the ARMA class. For * each method, the coefficients and forecasts for the last ten years, * 1860-1869, are compared. The method of moments and maximum likelihood * estimates produced similar results, but the least squares estimates were very * different from the other two.

* * @see Code * @see Output */ public class ARMAMaxLikelihoodEx1 { public static void main(String args[]) throws Exception { int backwardOrigin = 0, n_forecast = 10, n_series = 100; double[] sunspots = { 100.8, 81.6, 66.5, 34.8, 30.6, 7, 19.8, 92.5, 154.4, 125.9, 84.8, 68.1, 38.5, 22.8, 10.2, 24.1, 82.9, 132, 130.9, 118.1, 89.9, 66.6, 60, 46.9, 41, 21.3, 16, 6.4, 4.1, 6.8, 14.5, 34, 45, 43.1, 47.5, 42.2, 28.1, 10.1, 8.1, 2.5, 0, 1.4, 5, 12.2, 13.9, 35.4, 45.8, 41.1, 30.4, 23.9, 15.7, 6.6, 4, 1.8, 8.5, 16.6, 36.3, 49.7, 62.5, 67, 71, 47.8, 27.5, 8.5, 13.2, 56.9, 121.5, 138.3, 103.2, 85.8, 63.2, 36.8, 24.2, 10.7, 15, 40.1, 61.5, 98.5, 124.3, 95.9, 66.5, 64.5, 54.2, 39, 20.6, 6.7, 4.3, 22.8, 54.8, 93.8, 95.7, 77.2, 59.1, 44, 47, 30.5, 16.3, 7.3, 37.3, 73.9 }; double[] z; double arMM[], maMM[], constantMM; double arLS[], maLS[], constantLS; double ar[], ma[], constant; double forecastMM[][], forecastLS[][], forecast[][]; double deviations[]; double likelihood, var, varMM, varLS; double[] avgDev = {0.0, 0.0, 0.0}; ARMA armaMM, armaLS; ARMAMaxLikelihood maxArma; double[][] printOutput, printOutput2; String[] colLabels = { "Method of Moments", "Least Squares", "Maximum Likelihood" }; String[] colLabels1 = {"Least Squares", "Maximum Likelihood"}; String[] colLabels2 = { "Observed Sunspots", "Method of Moments", "Least Squares", "Maximum Likelihood" }; String[] colLabels3 = { "Lower Confidence Limit", "Forecast", "Upper Confidence Limit" }; PrintMatrixFormat pmf = new PrintMatrixFormat(); PrintMatrix pm = new PrintMatrix(); NumberFormat nf = NumberFormat.getNumberInstance(); pm.setColumnSpacing(3); nf.setMinimumFractionDigits(3); pmf.setNoRowLabels(); pmf.setNumberFormat(nf); pmf.setColumnLabels(colLabels); printOutput = new double[1][3]; printOutput2 = new double[n_forecast][4]; z = new double[n_series]; for (int i = 0; i < n_series; i++) { z[i] = sunspots[i]; } /* Method of Moments ARMA(2,1) Estimation */ armaMM = new ARMA(2, 1, z); armaMM.setMethod(ARMA.METHOD_OF_MOMENTS); armaMM.compute(); armaMM.setBackwardOrigin(backwardOrigin); arMM = armaMM.getAR(); maMM = armaMM.getMA(); constantMM = armaMM.getConstant(); forecastMM = armaMM.forecast(n_forecast); varMM = armaMM.getInnovationVariance(); /* Least Squares ARMA(2,1) Estimation */ armaLS = new ARMA(2, 1, z); armaLS.setMethod(ARMA.LEAST_SQUARES); armaLS.compute(); armaLS.setBackwardOrigin(backwardOrigin); arLS = armaLS.getAR(); maLS = armaLS.getMA(); constantLS = armaLS.getConstant(); varLS = armaLS.getInnovationVariance(); forecastLS = armaLS.forecast(n_forecast); /* Maximum Likelihood ARMA(2,1) Estimation */ maxArma = new ARMAMaxLikelihood(2, 1, z); maxArma.compute(); maxArma.setBackwardOrigin(backwardOrigin); ar = maxArma.getAR(); ma = maxArma.getMA(); constant = maxArma.getConstant(); likelihood = maxArma.getLikelihood(); var = maxArma.getInnovationVariance(); maxArma.setConfidence(0.9); forecast = maxArma.forecast(n_forecast); deviations = maxArma.getDeviations(); printOutput[0][0] = constantMM; printOutput[0][1] = constantLS; printOutput[0][2] = constant; pm.setTitle("ARMA(2,1) - Constant Term"); pm.print(pmf, printOutput); printOutput[0][0] = arMM[0]; printOutput[0][1] = arLS[0]; printOutput[0][2] = ar[0]; pm.setTitle("ARMA(2,1) - AR(1) Coefficient"); pm.print(pmf, printOutput); printOutput[0][0] = arMM[1]; printOutput[0][1] = arLS[1]; printOutput[0][2] = ar[1]; pm.setTitle("ARMA(2,1) - AR(2) Coefficient"); pm.print(pmf, printOutput); printOutput[0][0] = maMM[0]; printOutput[0][1] = maLS[0]; printOutput[0][2] = ma[0]; pm.setTitle("ARMA(2,1) - MA(1) Coefficient"); pm.print(pmf, printOutput); System.out.println("INNOVATION VARIANCE:"); System.out.println("Method of Moments " + varMM); System.out.println("Least Squares " + varLS); System.out.println("Maximum Likelihood " + var); System.out.println(""); for (int i = 0; i < n_forecast; i++) { printOutput2[i][0] = sunspots[90 + i]; printOutput2[i][1] = forecastMM[i][backwardOrigin]; printOutput2[i][2] = forecastLS[i][backwardOrigin]; printOutput2[i][3] = forecast[i][backwardOrigin]; } pm.setTitle("SUNSPOT FORECASTS FOR 1860-1869"); nf.setMaximumFractionDigits(0); pmf.setNumberFormat(nf); pmf.setColumnLabels(colLabels2); pmf.setFirstRowNumber(1860); pm.print(pmf, printOutput2); /* Get Confidence Interval Deviations */ printOutput2 = new double[n_forecast][3]; for (int i = 0; i < n_forecast; i++) { printOutput2[i][0] = Math.max(0, forecast[i][backwardOrigin] - deviations[i + backwardOrigin]); printOutput2[i][1] = forecast[i][backwardOrigin]; printOutput2[i][2] = forecast[i][backwardOrigin] + deviations[i]; } nf.setMaximumFractionDigits(0); pmf.setNumberFormat(nf); pmf.setColumnLabels(colLabels3); pmf.setFirstRowNumber(1860); pm.setTitle("SUNSPOT MAX. LIKELIHOOD 90% CONFIDENCE INTERVALS"); pm.print(pmf, printOutput2); } }