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 maximum likelihood estimates
are independently normally distributed with mean zero and variance 
. The maximum likelihood estimates 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. However, the least squares method does not converge using the default relative error convergence criteria. As a result, a warning is produced and the estimates are very different from the method of moments and maximum likelihood estimates.
using System;
using Imsl.Stat;
using Imsl.Math;
public class ARMAMaxLikelihoodEx1
{
public static void Main(String[] args)
{
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 = null;
double[] arMM, maMM;
double constantMM;
double[] arLS, maLS;
double constantLS;
double[] ar, ma;
double constant;
double[,] forecastMM, forecastLS, forecast;
double[] deviations;
double likelihood, var, varMM, varLS;
double[] avgDev = new double[]{0.0, 0.0, 0.0};
ARMA armaMM = null;
ARMA armaLS = null;
ARMAMaxLikelihood maxArma = null;
double[,] printOutput = null;
double[,] printOutput2 = null;
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();
pm.SetColumnSpacing(3);
pmf.SetNoRowLabels();
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.Method = Imsl.Stat.ARMA.ParamEstimation.MethodOfMoments;
armaMM.Compute();
armaMM.BackwardOrigin = backwardOrigin;
arMM = armaMM.GetAR();
maMM = armaMM.GetMA();
constantMM = armaMM.Constant;
forecastMM = armaMM.Forecast(n_forecast);
varMM = armaMM.InnovationVariance;
/* Least Squares ARMA(2,1) Estimation */
armaLS = new ARMA(2, 1, z);
armaLS.Method = Imsl.Stat.ARMA.ParamEstimation.LeastSquares;
armaLS.Compute();
armaLS.BackwardOrigin = backwardOrigin;
arLS = armaLS.GetAR();
maLS = armaLS.GetMA();
constantLS = armaLS.Constant;
varLS = armaLS.InnovationVariance;
forecastLS = armaLS.Forecast(n_forecast);
/* Maximum Likelihood ARMA(2,1) Estimation */
maxArma = new ARMAMaxLikelihood(2, 1, z);
maxArma.Compute();
maxArma.BackwardOrigin = backwardOrigin;
ar = maxArma.GetAR();
ma = maxArma.GetMA();
constant = maxArma.Constant;
likelihood = maxArma.Likelihood;
var = maxArma.InnovationVariance;
maxArma.Confidence = 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);
Console.Out.WriteLine("INNOVATION VARIANCE:");
Console.Out.WriteLine("Method of Moments " + varMM);
Console.Out.WriteLine("Least Squares " + varLS);
Console.Out.WriteLine("Maximum Likelihood " + var);
Console.Out.WriteLine("");
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");
pmf.SetColumnLabels(colLabels2);
pmf.FirstRowNumber = 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];
}
pmf.SetColumnLabels(colLabels3);
pmf.FirstRowNumber = 1860;
pm.SetTitle("SUNSPOT MAX. LIKELIHOOD 90% CONFIDENCE INTERVALS");
pm.Print(pmf, printOutput2);
}
}
ARMA(2,1) - Constant Term
Method of Moments Least Squares Maximum Likelihood
15.5439819435637 17.9316722242554 15.7580039183081
ARMA(2,1) - AR(1) Coefficient
Method of Moments Least Squares Maximum Likelihood
1.24425777984372 1.53128221858841 1.22506595122835
ARMA(2,1) - AR(2) Coefficient
Method of Moments Least Squares Maximum Likelihood
-0.575149766040151 -0.894433318954122 -0.56051392292258
ARMA(2,1) - MA(1) Coefficient
Method of Moments Least Squares Maximum Likelihood
-0.124089747872598 -0.132018760615022 -0.382891621968495
INNOVATION VARIANCE:
Method of Moments 287.242403738122
Least Squares 239.68797223859
Maximum Likelihood 214.508788425629
SUNSPOT FORECASTS FOR 1860-1869
Observed Sunspots Method of Moments Least Squares Maximum Likelihood
1860 95.7 87.6861734610048 98.0265069437306 87.9339016395945
1861 77.2 82.1446177467774 101.939296986813 82.0608538716595
1862 59.1 67.3203794962267 86.351331124405 66.9997857593059
1863 44 52.0624301952585 58.982026390731 51.8409090696475
1864 47 41.6037652345956 31.0142927589859 41.7122237493372
1865 30.5 37.3660959612162 12.6678176248844 37.8006776731389
1866 16.3 38.1086417046268 9.58945929412723 38.6860449014629
1867 7.3 41.469854513895 21.2853225650988 41.9631541830784
1868 37.3 45.2249946909406 41.9483762816881 45.4815685240649
1869 73.9 47.9641563097715 63.1281732161496 47.9549927562368
SUNSPOT MAX. LIKELIHOOD 90% CONFIDENCE INTERVALS
Lower Confidence Limit Forecast Upper Confidence Limit
1860 63.8431804305709 87.9339016395945 112.024622848618
1861 36.4438944595113 82.0608538716595 127.677813283808
1862 10.134623371631 66.9997857593059 123.864948146981
1863 0 51.8409090696475 112.081240519376
1864 0 41.7122237493372 102.18743849516
1865 0 37.8006776731389 98.4521257842839
1866 0 38.6860449014629 99.9504700070846
1867 0 41.9631541830784 103.74778523107
1868 0 45.4815685240649 107.46464011018
1869 0 47.9549927562368 109.958371498445
Imsl.Stat.ARMA: Relative function convergence - Both the scaled actual and
predicted reductions in the function are less than or equal to the relative
function convergence tolerance "ConvergenceTolerance" = 5.00000041370186E-10.
Imsl.Stat.ARMA: Least squares estimation of the parameters has failed to
converge. Increase "length" and/or "tolerance" and/or "convergence_tolerance".
The estimates of the parameters at the last iteration may be used as new
starting values.
Link to C# source.