Example 2 (Forecasting): ARAutoUnivariate

Using the Canadian Lynx data included in TIMSAC-78, ARAutoUnivariate is used to find the minimum AIC autoregressive model using a maximum number of lags of maxlag =20.

This example compares the three different methods for estimating the autoregressive coefficients, and it illustrates the relationship between these estimates and those calculated within the TIMSAC UNIMAR code. As illustrated, the UNIMAR code estimates the coefficients and innovation variance using only the last N-maxlag values in the time series. The other estimation methods use all N-k values, where k is the number of lags with minimum AIC selected by ARAutoUnivariate .

This example also illustrates how to generate forecasts for the observed series values and beyond by setting the backward orgin for the forecasts.

using System;
using Imsl.Stat;
using Imsl.Math;
using PrintMatrix = Imsl.Math.PrintMatrix;	
public class ARAutoUnivariateEx2
{
	public static void  Main(String[] args)
	{
		/* THE CANDIAN LYNX DATA AS USED IN TIMSAC 1821-1934 */
		double[] y = new double[]{
                   0.24300e01, 0.25060e01, 0.27670e01, 0.29400e01, 0.31690e01,
                   0.34500e01, 0.35940e01, 0.37740e01, 0.36950e01, 0.34110e01, 
                   0.27180e01, 0.19910e01, 0.22650e01, 0.24460e01, 0.26120e01, 
                   0.33590e01, 0.34290e01, 0.35330e01, 0.32610e01, 0.26120e01, 
                   0.21790e01, 0.16530e01, 0.18320e01, 0.23280e01,  0.27370e01, 
                   0.30140e01, 0.33280e01, 0.34040e01, 0.29810e01, 0.25570e01,  
                   0.25760e01, 0.23520e01, 0.25560e01, 0.28640e01, 0.32140e01, 
                   0.34350e01,  0.34580e01, 0.33260e01, 0.28350e01, 0.24760e01, 
                   0.23730e01, 0.23890e01,  0.27420e01, 0.32100e01, 0.35200e01, 
                   0.38280e01, 0.36280e01, 0.28370e01,  0.24060e01, 0.26750e01, 
                   0.25540e01, 0.28940e01, 0.32020e01, 0.32240e01,  0.33520e01, 
                   0.31540e01, 0.28780e01, 0.24760e01, 0.23030e01, 0.23600e01,  
                   0.26710e01, 0.28670e01, 0.33100e01, 0.34490e01, 0.36460e01, 
                   0.34000e01,  0.25900e01, 0.18630e01, 0.15810e01, 0.16900e01, 
                   0.17710e01, 0.22740e01,  0.25760e01, 0.31110e01, 0.36050e01, 
                   0.35430e01, 0.27690e01, 0.20210e01,  0.21850e01, 0.25880e01, 
                   0.28800e01, 0.31150e01, 0.35400e01, 0.38450e01,  0.38000e01, 
                   0.35790e01, 0.32640e01, 0.25380e01, 0.25820e01, 0.29070e01,  
                   0.31420e01, 0.34330e01, 0.35800e01, 0.34900e01, 0.34750e01, 
                   0.35790e01,  0.28290e01, 0.19090e01, 0.19030e01, 0.20330e01, 
                   0.23600e01, 0.26010e01, 0.30540e01, 0.33860e01, 0.35530e01, 
                   0.34680e01, 0.31870e01, 0.27230e01, 0.26860e01, 0.28210e01, 
                   0.30000e01, 0.32010e01, 0.34240e01, 0.35310e01};
		double[][] printOutput = null;
		double[] timsacAR, mmAR, mleAR, lsAR;
		double[] forecasts, residuals;
		double timsacConstant, mmConstant, mleConstant, lsConstant;
		double timsacVar, timsacEquivalentVar, mmVar, mleVar, lsVar;
		int maxlag = 20;
		String[] colLabels = new String[]{"TIMSAC", 
                                          "Method of Moments", 
                                          "Least Squares", 
                                          "Maximum Likelihood"};
		String[] colLabels2 = new String[]{"Observed", "Forecast", "Residual"};
		PrintMatrixFormat pmf = new PrintMatrixFormat();
		PrintMatrix pm = new PrintMatrix();
		pmf.SetColumnLabels(colLabels);
        pmf.NumberFormat = "0.0000";
		Console.Out.WriteLine
            ("Automatic Selection of Minimum AIC AR Model");
		Console.Out.WriteLine("");
		
		ARAutoUnivariate autoAR = new ARAutoUnivariate(maxlag, y);
		autoAR.Compute();
		int orderSelected = autoAR.Order;
		Console.Out.WriteLine("Minimum AIC Selected=" + autoAR.AIC 
            + " with an optimum lag of k= " + autoAR.Order);
		Console.Out.WriteLine("");
		
		timsacAR = autoAR.GetTimsacAR();
		timsacConstant = autoAR.TimsacConstant;
		timsacVar = autoAR.TimsacVariance;
		lsAR = autoAR.GetAR();
		lsConstant = autoAR.Constant;
		lsVar = autoAR.InnovationVariance;
		
		autoAR.EstimationMethod = 
            Imsl.Stat.ARAutoUnivariate.ParamEstimation.MethodOfMoments;
		autoAR.Compute();
		mmAR = autoAR.GetAR();
		mmConstant = autoAR.Constant;
		mmVar = autoAR.InnovationVariance;
		
		autoAR.EstimationMethod = 
            Imsl.Stat.ARAutoUnivariate.ParamEstimation.MaximumLikelihood;
		autoAR.Compute();
		mleAR = autoAR.GetAR();
		mleConstant = autoAR.Constant;
		mleVar = autoAR.InnovationVariance;
		
		printOutput = new double[orderSelected + 1][];
		for (int i = 0; i < orderSelected + 1; i++)
		{
			printOutput[i] = new double[4];
		}
		printOutput[0][0] = timsacConstant;
		for (int i = 0; i < orderSelected; i++)
			printOutput[i + 1][0] = timsacAR[i];
		printOutput[0][1] = mmConstant;
		for (int i = 0; i < orderSelected; i++)
			printOutput[i + 1][1] = mmAR[i];
		printOutput[0][2] = lsConstant;
		for (int i = 0; i < orderSelected; i++)
			printOutput[i + 1][2] = lsAR[i];
		printOutput[0][3] = mleConstant;
		for (int i = 0; i < orderSelected; i++)
			printOutput[i + 1][3] = mleAR[i];
		pm.SetTitle("Comparison of AR Estimates");
		pm.Print(pmf, printOutput);
		
		/* calculation of equivalent innovation variance using TIMSAC 
		coefficients.  The Timsac innovation variance is calculated using
		only N-maxlag observations in the series.  The following code 
		calculates the innovation variance using N-k observations in the 
		series with the Timsac coefficient.  This illustrates that the
		least squares Timsac coefficients will not have the least value for
		the sum of squared residuals, which is calculated using all N-k
		observations. */
		ARMA armaLS = new ARMA(orderSelected, 0, y);
		armaLS.SetARMAInfo(autoAR.TimsacConstant, autoAR.GetTimsacAR(), 
            new double[0], autoAR.TimsacVariance);
		armaLS.BackwardOrigin = y.Length - orderSelected;
		forecasts = armaLS.GetForecast(1);
		double sumResiduals = 0.0;
		for (int i = 0; i < y.Length - orderSelected; i++)
		{
			sumResiduals += (y[i + orderSelected] - forecasts[i]) * 
                (y[i + orderSelected] - forecasts[i]);
		}
		timsacEquivalentVar = sumResiduals / (y.Length - orderSelected - 1);
		printOutput = new double[1][];
		for (int i2 = 0; i2 < 1; i2++)
		{
			printOutput[i2] = new double[4];
		}
		printOutput[0][0] = timsacEquivalentVar;
		/* the method of moments variance */
		printOutput[0][1] = mmVar;
		/* the least squares variance */
		printOutput[0][2] = lsVar;
		/* the maximum likelihood estimate of the variance */
		printOutput[0][3] = mleVar;
		pm.SetTitle("Comparison of Equivalent Innovation Variances");
		pm.Print(pmf, printOutput);
		
		/* FORECASTING - An example of forecasting using the maximum
		* likelihood estimates for the minimum AIC AR model.  In this example,
		* forecasts are returned for the last 10 values in the series followed
		* by the forecasts for the next 5 values.
		*/
		autoAR.BackwardOrigin = 10;
		forecasts = autoAR.GetForecast(15);
		residuals = autoAR.GetResiduals();
		printOutput = new double[15][];
		for (int i3 = 0; i3 < 15; i3++)
		{
			printOutput[i3] = new double[3];
		}
		for (int i = 0; i < 10; i++)
		{
			printOutput[i][0] = y[y.Length - 10 + i];
			printOutput[i][1] = forecasts[i];
			printOutput[i][2] = residuals[i];
		}
		for (int i = 10; i < 15; i++)
		{
			printOutput[i][0] = Double.NaN;
			printOutput[i][1] = forecasts[i];
			printOutput[i][2] = Double.NaN;
		}
		pmf.FirstRowNumber = 105;
		pmf.SetColumnLabels(colLabels2);
		pm.SetTitle("Maximum Likelihood Forecasts of Last 10 Values");
		pm.Print(pmf, printOutput);
	}
}

Output

Automatic Selection of Minimum AIC AR Model

Minimum AIC Selected=-296.130132635624 with an optimum lag of k= 11

                    Comparison of AR Estimates
    TIMSAC   Method of Moments  Least Squares  Maximum Likelihood  
 0   1.0427        1.1679           1.1144           1.1189        
 1   1.1813        1.1381           1.1481           1.1664        
 2  -0.5516       -0.5061          -0.5331          -0.5419        
 3   0.2314        0.2098           0.2757           0.2624        
 4  -0.1780       -0.2672          -0.3263          -0.3052        
 5   0.0199        0.1112           0.1685           0.1519        
 6  -0.0626       -0.1246          -0.1643          -0.1460        
 7   0.0286        0.0693           0.0728           0.0581        
 8  -0.0507       -0.0419          -0.0305          -0.0310        
 9   0.1999        0.1366           0.1509           0.1380        
10   0.1618        0.1828           0.1935           0.1995        
11  -0.3391       -0.3101          -0.3414          -0.3376        

          Comparison of Equivalent Innovation Variances
   TIMSAC  Method of Moments  Least Squares  Maximum Likelihood  
0  0.0377       0.0427           0.0369            0.0362        

Maximum Likelihood Forecasts of Last 10 Values
     Observed  Forecast  Residual  
105   3.5530    3.4388    0.1142   
106   3.4680    3.4801   -0.0121   
107   3.1870    2.9243    0.2627   
108   2.7230    2.7026    0.0204   
109   2.6860    2.5558    0.1302   
110   2.8210    2.7852    0.0358   
111   3.0000    2.9493    0.0507   
112   3.2010    3.1861    0.0149   
113   3.4240    3.3855    0.0385   
114   3.5310    3.5272    0.0038   
115   NaN       3.4465    NaN      
116   NaN       3.1947    NaN      
117   NaN       2.8289    NaN      
118   NaN       2.4917    NaN      
119   NaN       2.4142    NaN      


Link to C# source.