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, but the least squares estimates were very different from the other two.
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;
Imsl.WarningObject w = Imsl.Warning.WarningObject;
Imsl.Warning.WarningObject = null;
armaLS.Compute();
Imsl.Warning.WarningObject = w;
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.7580039181538
ARMA(2,1) - AR(1) Coefficient
Method of Moments Least Squares Maximum Likelihood
1.24425777984372 1.53128221858841 1.22506595122502
ARMA(2,1) - AR(2) Coefficient
Method of Moments Least Squares Maximum Likelihood
-0.575149766040151 -0.894433318954122 -0.560513922915966
ARMA(2,1) - MA(1) Coefficient
Method of Moments Least Squares Maximum Likelihood
-0.124089747872598 -0.132018760615022 -0.382891621972049
INNOVATION VARIANCE:
Method of Moments 287.242403738122
Least Squares 239.68797223859
Maximum Likelihood 214.508788425648
SUNSPOT FORECASTS FOR 1860-1869
Observed Sunspots Method of Moments Least Squares Maximum Likelihood
1860 95.7 86.0415456007167 97.7310653809498 85.3832083889174
1861 77.2 80.0982767366663 101.486892575095 78.9360864182773
1862 59.1 65.7201111032504 85.9228250706849 64.6014386269177
1863 44 51.2482363506636 58.7305082695945 50.6542513222013
1864 47 41.5110922009792 31.0124176242576 41.602796706989
1865 30.5 37.7190702173874 12.8899124522828 38.3317605186442
1866 16.3 38.6011335425551 9.93122633714127 39.3979917935465
1867 7.3 41.8796282538917 21.6100153473733 42.5376567504899
1868 37.3 45.451602289546 42.1398847350968 45.7863159128644
1869 73.9 48.0104333066195 63.1310106626735 48.0063117181664
SUNSPOT MAX. LIKELIHOOD 90% CONFIDENCE INTERVALS
Lower Confidence Limit Forecast Upper Confidence Limit
1860 61.2924871798927 85.3832083889174 109.473929597942
1861 33.3191270061225 78.9360864182773 124.553045830432
1862 7.73627623921441 64.6014386269177 121.466601014621
1863 0 50.6542513222013 110.894582772018
1864 0 41.602796706989 102.078011452932
1865 0 38.3317605186442 98.9832086298742
1866 0 39.3979917935465 100.662416899195
1867 0 42.5376567504899 104.322287798477
1868 0 45.7863159128644 107.769387498973
1869 0 48.0063117181664 110.009690460371
Link to C# source.