Products & Services Solutions Academia Support User Community Company

Learn more about Econometrics Toolbox   

Model Identification

On this page…

Autocorrelation and Partial Autocorrelation

Unit Root Tests

Misspecification Tests

Akaike and Bayesian Information

Autocorrelation and Partial Autocorrelation

See Example: Using the Default Model for information about using the autocorrelation (autocorr) and partial autocorrelation (parcorr) functions as qualitative guides in the process of model selection and assessment. This example also introduces the following functions:

Back to Top

Unit Root Tests

What Is a Unit Root Test?

A unit root process is a data-generating process whose first difference is stationary. In other words, a unit root process yt has the form

yt = yt–1 + stationary process.

A unit root test attempts to determine whether a given time series is consistent with a unit root process.

The next section gives more details of unit root processes, and suggests why it is important to detect them.

Modeling Unit Root Processes

There are two basic models for economic data with linear growth characteristics:

Here Δ is the differencing operator, Δyt = yt – yt–1.

The processes are indistinguishable for finite data. In other words, there are both a TSP and a DSP that fit a finite data set arbitrarily well. However, the processes are distinguishable when restricted to a particular subclass of data-generating processes, such as AR(p) processes. After fitting a model to data, a unit root test checks if the AR(1) coefficient is 1.

There are two main reasons to distinguish between these types of processes:

Forecasting.   A TSP and a DSP produce different forecasts. Basically, shocks to a TSP return to the trend line c + dt as time increases. In contrast, shocks to a DSP might be persistent over time.

For example, consider the simple trend-stationary model

y1(t) = 0.9y1(t–1) + 0.02t + ε1(t)

and the difference-stationary model

y2(t) = 0.2 + y2(t–1) + ε2(t).

In these models, ε1(t) and ε2(t) are independent innovation processes. For this example, the innovations are independent and distributed N(0,1).

Both processes grow at rate 0.2t. To calculate the growth rate for the TSP, which has a linear term 0.02t, set ε1(t) = 0. Then solve the model y1(t) = a + ct for a and c:

a + ct = 0.9(a + c(t–1)) + 0.02t.

The solution is a = –1.8, c = 0.2.

A plot for t = 1:1000 shows the TSP stays very close to the trend line, while the DSP has persistent deviations away from the trend line.

Forecasts based on the two series are different. To see this difference, plot the predicted behavior of the two series using the vgxpred function. The following plot shows the last 100 data points in the two series and predictions of the next 100 points, including confidence bounds.

 Code for Generating the Figure

The TSP has confidence intervals that do not grow with time. Whereas, the DSP has confidence intervals that grow. Furthermore, the TSP goes to the trend line quickly, while the DSP does not tend towards the trend line y = 0.2t asymptotically.

Spurious Regression.   The presence of unit roots can lead to false inferences in regressions between time series.

Suppose xt and yt are unit root processes with independent increments, such as random walks with drift

xt = c1 + xt–1 + ε1(t)
yt = c2 + yt–1 + ε2(t),

where εi(t) are independent innovations processes. Regressing y on x results, in general, in a nonzero regression coefficient, and significant coefficient of determination R2. This result holds despite xt and yt being independent random walks.

If both processes have trends (ci ≠ 0), there is a correlation between x and y because of their linear trends. However, even if the ci = 0, the presence of unit roots in the xt and yt processes yields correlation. For more information on spurious regression, see Granger and Newbold [27].

Available Tests

There are four Econometrics Toolbox tests for unit roots. These functions test for the existence of a single unit root; when there are two or more unit roots, their results might not be valid.

Dickey-Fuller and Phillips-Perron Tests.   adftest performs the augmented Dickey-Fuller test. pptest performs the Phillips-Perron test. These two classes of tests have a null hypothesis of a unit root process of the form

yt = yt–1 + c + dt + et,

which the functions test against an alternative model

yt = ayt–1 + c + dt + et,

where a < 1. The null and alternative models for a Dickey-Fuller test are like those for a Phillips-Perron test. The difference is adftest extends the model with extra parameters accounting for serial correlation among the residues:

yt = c + dt + ayt – 1 + b1Δyt – 1 + b2Δyt – 2 +...+ bpΔyt – p + e(t),

where

Phillips-Perron adjusts the test statistics to account for serial correlation.

There are three variants of both adftest and pptest, corresponding to the following values of the 'model' parameter:

For information on how to choose the appropriate value of 'model', see Choose Models to Test.

KPSS Test.   The KPSS test, kpsstest, is an inverse of the Phillips-Perron test: it reverses the null and alternative hypotheses. The KPSS test uses the model:

yt = ct + dt + ut, with
ct = ct–1 + vt.

Here ut is a stationary process, and vt is an i.i.d. process with mean 0 and variance σ2. The null hypothesis is that σ2 = 0, so that the random walk term ct becomes a constant intercept. The alternative is σ2 > 0, which introduces the unit root in the random walk.

Variance Ratio Test.   The variance ratio test, vratiotest, is based on the fact that the variance of a random walk increases linearly with time. vratiotest can also take into account heteroskedasticity, where the variance increases at a variable rate with time. The test has a null hypotheses of a random walk:

Δyt = et.

Testing for Unit Roots

Transform Data.   Transform your time series to be approximately linear before testing for a unit root. If a series has exponential growth, take its logarithm. For example, GDP and consumer prices typically have exponential growth, so test their logarithms for unit roots. There is more discussion of this topic in Data Preprocessing.

If you want to transform your data to be stationary instead of approximately linear, unit root tests can help you determine whether to difference your data, or to subtract a linear trend. For a discussion of this topic, see Modeling Unit Root Processes.

Choose Models to Test.  

Determine Appropriate Lags.   Setting appropriate lags depends on the test you use:

Run Tests.   Run multiple tests simultaneously by entering a vector of parameters for lags, alpha, model, or test. All vector parameters must have the same length. The test expands any scalar parameter to the length of a vector parameter. For an example using this technique, see Example: Testing Wage Data for a Unit Root.

Examples of Unit Root Tests

Example: Testing Simulated Data for a Unit Root.   This example generates trend-stationary, difference-stationary, stationary (AR(1)), and heteroskedastic random walk time series using simulated data. The example shows how to test the resulting series. It also shows that the tests give the expected results.

  1. Generate four time series:

    T = 1e3; % Sample size
strm = RandStream('mt19937ar','Seed',142857); % reproducible
t = (1:T)'; % time multiple

y1 = randn(strm,T,1) + .2*t; % trend stationary
y2 = randn(strm,T,1) + .2;
y2 = cumsum(y2); % difference stationary
y3 = randn(strm,T,1) + .2;
for ii = 2:T
    y3(ii) = y3(ii) + .99*y3(ii-1); % AR(1)
end
sigma = (sin(t/200) + 1.5)/2; % std deviation
y4 = randn(strm,T,1).*sigma + .2;
y4 = cumsum(y4); % heteroskedastic

  2. Plot the first 100 points in each series:

    y = [y1 y2 y3 y4];
plot1 = plot(y(1:100,:));
set(plot1(1),'LineWidth',2)
set(plot1(3),'LineStyle',':','LineWidth',2)
set(plot1(4),'LineStyle',':','LineWidth',2)
legend('trend stationary','difference stationary','AR(1)',...
   'Heteroskedastic','location','northwest')

    It is difficult to distinguish which series comes from which model simply by looking at these initial segments of the series.

  3. Plot the entire data set:

    plot2 = plot(y);
set(plot2(1),'LineWidth',2)
set(plot2(3),'LineStyle',':','LineWidth',2)
set(plot2(4),'LineStyle',':','LineWidth',2)
legend('trend stationary','difference stationary','AR(1)',...
   'Heteroskedastic','location','northwest')

    The differences between the series are clearer here:

    • The trend stationary series has little deviation from its mean trend.

    • The difference stationary and Heteroskedastic series have persistent deviations away from the trend line.

    • The AR(1) series exhibits long-run stationary behavior; the others grow linearly.

    • The difference stationary and Heteroskedastic series appear similar. But notice that the Heteroskedastic series has much more local variability near time 300, and much less near time 900. The model variance is maximal when sin(t/200) = 1, at time 200*π/2 ≈ 314. The model variance is minimal when sin(t/200) = –1, at time 200*3*π/2 ≈ 942. So the visual variability matches the model.

  4. Perform the following tests (the number of lags is set to 2 for demonstration purposes):

    • Test the three growing series:

      • Null: yt = yt – 1 + c + b1Δyt – 1 + b2Δyt–2 + e(t)

      • Alternative: yt = ayt – 1 + c + dt + b1Δyt – 1 + b2Δyt–2 + e(t)

    • Test the stationary AR(1) series:

      • Null: yt = yt – 1 + b1Δyt – 1 + b2Δyt–2 + e(t)

      • Alternative: yt = ayt – 1 + b1Δyt – 1 + b2Δyt–2 + e(t)

    • Test all the series:

      • Null: Var(Δyt) = constant

      • Alternative: Var(Δyt) ≠ constant

    Test Results

    Tested ModelData Generation Process
    Trend stationary y1Difference stationary y2Stationary AR(1) y3Heteroskedastic Random Walk y4
    'TS'

    adftest(y1, 'model','ts', 'lags',2)

    1

    adftest(y2, 'model','ts', 'lags',2)

    0

    		 

    adftest(y4, 'model','ts', 'lags',2)

    0

    'AR'  

    adftest(y3, 'model','ar', 'lags',2)

    0

    		 
    Stationary (KPSS)

    kpsstest(y1, 'lags',32, 'trend',true)

    0

    kpsstest(y2, 'lags',32, 'trend',true)

    1

    kpsstest(y3, 'lags',32)

    1

    kpsstest(y4, 'lags',32, 'trend',true)

    1

    Random Walk

    vratiotest(y1)

    1

    vratiotest(y2, 'IID',true)

    0

    vratiotest(y3, 'IID',true)

    0

    vratiotest(y4)

    0

    vratiotest (y4, 'IID',true)

    0 ×

    In each cell of the table:

    • The first line represents the command.

    • The next line represents the result of running the command.

    • The final symbol indicates whether the result matches the expected result based on the data generation process: for a match, × for failure to match.

    • You can substitute pptest for adftest.

    There is one incorrect inference. In the lower right corner of the table, vratiotest does not reject the hypothesis that the heteroskedastic process is an IID random walk. This inference depends on the pseudorandom numbers. Run the data generation procedure with Seed equal to 857142 and the inference changes.

Example: Testing Wage Data for a Unit Root.   This example tests whether there is a unit root in wages in the manufacturing sector. The data for the years 1900–1970 is in the Nelson-Plosser data set [45].

  1. Load the Nelson-Plosser data. Extract the ninth column (manufacturing wages):

    load NelsonPlosser
wages = NPData(:,9);

  2. Trim the NaN data from the series and the associated dates. (This step is optional, since the test ignores NaN values.)

    wdates = dates(isfinite(wages));
wages = wages(isfinite(wages));

  3. Plot the data to look for trends:

    plot(wdates,wages)
title('Wages')

  4. The plot suggests exponential growth. Use the log of wages to linearize the data set:

    lwages = log(wages);
plot(wdates,lwages)
title('Log Wages')

  5. The data appear to have a linear trend. Test the hypothesis that this process is a unit root process with a trend (difference stationary), against the alternative that there is no unit root (trend stationary). Following Determine Appropriate Lags, set 'lags' = [2:4]:

    [h,pValue] = pptest(lwages,'lags',[2:4],'model','TS')

h =
     0     0     0

pValue =
    0.6303    0.6008    0.5913

    pptest fails to reject the hypothesis of a unit root at all lags.

  6. Check the inference using kpsstest. Following Determine Appropriate Lags, set 'lags' to [7:2:11]:

    [h,pValue] = kpsstest(lwages,'lags',[7:2:11])
Warning: Test statistic #1 below tabulated critical values:
maximum p-value = 0.100 reported. 
> In kpsstest>getStat at 608
  In kpsstest at 292

h =
     0     0     0

pValue =
    0.1000    0.1000    0.1000

    kpsstest fails to reject the hypothesis that the data is trend-stationary. If the result would have been [1 1 1], the two inferences would provide consistent evidence of a unit root. It remains unclear whether the data has a unit root. This is a typical result of tests on many macroeconomic series.

    The warning that the test statistic "...is below tabulated critical values" does not indicate a problem. kpsstest has a limited set of calculated critical values. When the calculated test statistic is outside this range, the test reports the p-value at the endpoint. So, in this case, pValue reflects the closest tabulated value. When a test statistic lies inside the span of tabulated values, kpsstest reports linearly interpolated values.

Example: Testing Stock Data for a Random Walk.   This example uses market data for daily returns of stocks and cash (money market) from the period January 1, 2000 to November 7, 2005. The question is whether a random walk is a good model for the data.

  1. Load the data.

    load CAPMuniverse

  2. Extract two series to test. The first column of data is the daily return of a technology stock. The last (14th) column is the daily return for cash (the daily money market rate).

    tech1 = Data(:,1);
money = Data(:,14);

    The returns are the logs of the ratios of values at the end of a day over the values at the beginning of the day.

  3. Convert the data to prices (values) instead of returns. vratiotest takes prices as inputs, as opposed to returns.

    tech1 = cumsum(tech1);
money = cumsum(money);

  4. Plot the data to see whether they appear to be stationary:

    subplot(2,1,1)
plot(Dates,tech1);
title('Log(relative stock value)')
datetick('x')
hold('on')
subplot(2,1,2);
plot(Dates,money)
title('Log(accumulated cash)')
datetick('x')
hold('off')

    Cash has a small variability, and appears to have long-term trends. The stock series has a good deal of variability, and no definite trend, though it appears to increase towards the end.

  5. Test whether the stock series matches a random walk:

    [h,pValue,stat,cValue,ratio] = vratiotest(tech1)

h =
     0

pValue =
    0.1646

stat =
   -1.3899

cValue =
    1.9600

ratio =
    0.9436

    vratiotest does not reject the hypothesis that a random walk is a reasonable model for the stock series.

  6. Test whether an i.i.d. random walk is a reasonable model for the stock series:

    [h,pValue,stat,cValue,ratio] = vratiotest(tech1,'IID',true)

h =
     1

pValue =
    0.0304

stat =
   -2.1642

cValue =
    1.9600

ratio =
    0.9436

    vratiotest rejects the hypothesis that an i.i.d. random walk is a reasonable model for the tech1 stock series at the 5% level. Thus, vratiotest indicates that the most appropriate model of the tech1 series is a heteroskedastic random walk.

  7. Test whether the cash series matches a random walk:

    [h,pValue,stat,cValue,ratio] = vratiotest(money)

h =
     1

pValue =
  4.6093e-145

stat =
   25.6466

cValue =
    1.9600

ratio =
    2.0006

    vratiotest emphatically rejects the hypothesis that a random walk is a reasonable model for the cash series (pValue = 4.6093e-145). Removing a linear trend from the cash series does not affect the resulting statistics.

Back to Top

Misspecification Tests

The Classical Misspecification Tests

Econometrics Toolbox includes three classical frameworks for model misspecification testing:

These functions all estimate the difference in maximum loglikelihood between a restricted model and an unrestricted model. They test whether unrestricted model parameters are statistically significant.

Some examples of model restrictions that can be tested in these frameworks are:

The setting for all of these tests is a vector of model parameters θ, data d, and a loglikelihood function L(θ|d). When data represent independent observations, L is a sum of individual loglikelihoods over all data points. The unrestricted optimizer maximizes L. The optimizer subject to restrictions r(θ) = 0 occurs, in general, at some . Therefore,

Each test assesses the difference in loglikelihood in a different way:

lratiotest evaluates the difference in loglikelihoods directly. lmtest and waldtest do so indirectly, with the idea that insignificant changes in the evaluated quantities are equivalent to insignificant changes in the parameters. This equivalence depends on the curvature of the loglikelihood surface in the neighborhood of the maximum likelihood estimator. As a result, the Wald and Lagrange multiplier tests include an estimate of parameter covariance in the formulation of the test statistic.

Each test requires different input data. The following table shows the input requirements for the tests, and some MATLAB functions that can provide the inputs.

Types of Hypothesis Tests

TestRequired InputsHelper Functions for the Inputs
lratiotestUnrestricted maximum loglikelihood (uLL)

garchfit, vgxvarx, Optimization Toolbox fmincon, fminunc, Statistics Toolbox mle, betalike, explike, and other functions named *like (see Negative Log-Likelihood Functions in the Statistics Toolbox user guide documentation).

Restricted maximum loglikelihood (rLL)

garchfit, vgxvarx, Optimization Toolbox fmincon

Degrees of freedom

Number of restrictions

lmtestScore at restricted maximum

Optimization Toolbox fmincon

Estimated unrestricted covariance at restricted maximum

You must compute the Jacobian analytically

Degrees of freedom

Number of restrictions

waldtestRestriction function r at unrestricted maximum

Optimization Toolbox fmincon, fminunc can find the maximum point

R = Jacobian of r at unrestricted maximum

Optimization Toolbox fmincon, fminunc can find the maximum point; you must compute the Jacobian analytically

Estimated unrestricted covariance at unrestricted maximum

garchfit, vgxvarx, Statistics Toolbox mlecov

Covariance Estimators

lmtest and waldtest require an estimate of the parameter covariance. Parameter covariance is related to the curvature of the loglikelihood function. This curvature indirectly indicates parameter significance used by the two tests.

There are several common methods for estimating the parameter covariance, including:

For lmtest, evaluate the estimator at the restricted maximum point . For waldtest, evaluate the estimator at the unrestricted maximum point .

garchfit and vgxvarx return an estimated covariance matrix, along with other parameter estimates. garchfit returns an estimator based on OPG, and vgxvarx returns an estimator based on the expected Hessian.

Example: Specifying Distribution Parameters

The following reproduces the example in Greene [28], Sec. 16.6.4. The example fits a Gamma distribution to the simulated education and income data in IncomeData.mat:

load IncomeData
x = education;
y = income;

The unrestricted model has a likelihood l given by the probability density function for the Gamma distribution:

where

This yields the loglikelihood

(3-1)

Gamma distributions are sums of ρ exponential distributions, so admit natural restrictions on the value of ρ. Greene's restricted model is ρ = 1, a simple exponential distribution.

lratiotest.   To perform this test, you need estimates of both the unrestricted and restricted maximum likelihood. To find the unrestricted maximum likelihood estimates of β and ρ, you can use fminunc:

% Negative loglikelihood function with parameters p = [beta,rho]:
nLLF = @(p)sum(p(2)*(log(p(1) + x)) + gammaln(p(2))...
    -(p(2) - 1)*log(y) + y./(p(1) + x));

options = optimset('LargeScale','off'); % 'on' with gradients
% Unrestricted estimate: General gamma model
[umle,unLL] = fminunc(nLLF,[.1 .1],options)

umle =
   -4.7186    3.1509

unLL =
   82.9160

ubeta = umle(1);    % Unrestricted beta estimate
urho = umle(2);     % Unrestricted rho estimate
uLL = -unLL;        % Unrestricted loglikelihood

lratiotest also requires the restricted maximum likelihood:

% Restricted estimate with rho = p(2) = 1
rnLLF = @(p)nLLF([p 1]);
[rmle,rnLL] = fminunc(rnLLF,.1,options)

rmle =
   15.6027

rnLL =
   88.4363

rbeta = rmle(1);    % Restricted beta estimate
rrho = 1;           % Restricted rho
rLL = -rnLL;        % Restricted loglikelihood

To evaluate the likelihood ratio, be careful to use the negative of unLL and rnLL. There is one restriction, sodof = 1:

dof = 1;
[LRh,LRp,LRstat,cV] = lratiotest(uLL,rLL,dof)

LRh =
     1

LRp =
  8.9146e-004

LRstat =
   11.0404

cV =
    3.8415

LRh = 1 means lratiotest rejects the restricted hypothesis. This is a strong rejection: the p-value is less than 1/1000.

lmtest.   lmtest requires the score, and an estimate of the covariance matrix, evaluated at the restricted parameter estimates.

By analytic calculation, the gradient of the loglikelihood function Equation 3-1 has components:

Here Ψ(ρ) is the digamma (psi) function, the derivative of the logarithm of Γ(ρ), and n is the number of terms in the sum, which is the number of data points.

The OPG estimator of the covariance is (G'G)–1, where G is the matrix of gradients. For more information, see Covariance Estimators. Since the score is evaluated as the sum of these derivatives, the simplest way to estimate the covariance for lmtest is to use the OPG estimator.

Calculate the gradient of L at the restricted maximum β = rbeta, ρ = 1, and calculate the OPG estimator, as follows:

RCG = [-rrho./(rbeta+x)+y.*(rbeta+x).^(-2),...
     -log(rbeta+x)-psi(rrho)+log(y)];
Rscore = sum(RCG)';
REstCov1 = inv(RCG'*RCG);

There is one restriction, so dof = 1.

[LMh,LMp,LMstat,LMcv] = lmtest(Rscore,REstCov1,dof)

LMh =
     1

LMp =
  7.4744e-005

LMstat =
   15.6868

LMcv =
    3.8415

LMh = 1 means lmtest rejects the restricted hypothesis. This is a strong rejection: the p-value is less than 1/10000.

waldtest.   waldtest requires the value and gradient of the restriction function at the unrestricted maximum. The restriction function is r(β, ρ) = ρ – 1 = 0. This has gradient [0 1]. In other words, at the unrestricted maximum  = (β, ρ),

r = 3.1509 – 1
R = [0 1].

R = [0 1];r = urho - 1;n = size(x,1);

The Wald test also requires an estimated covariance. This time, use the sum of observed Hessians,

which you find through analytic calculation:

(3-2)

The estimated covariance is given by –H–1, where H is the observed Hessian; see Covariance Estimators. Calculate the Hessian and covariance matrices as follows:

UDPsi = psi(1,urho);
UH = [sum(urho./(ubeta+x).^2)-2*sum(y./(ubeta+x).^3),...
     -sum(1./(ubeta+x)); -sum(1./(ubeta+x)),-n*UDPsi];
UEstCov2 = -inv(UH);

The Wald test is

[Wh,Wp,Wstat,Wcv] = waldtest(r,R,UEstCov2)

Wh =
     1

Wp =
    0.0068

Wstat =
    7.3337

Wcv =
    3.8415

Wh = 1 means waldtest rejects the restricted hypothesis.

Summary of Calculations for the Tests.   To compare these tests with Greene [28], this section collects the inputs for the three tests: the unrestricted and restricted values of the loglikelihoods, β, ρ, gradients, and covariance estimates. The results appear in the following tables.

QuantityVariableUnrestricted EstimateRestricted Estimate
βubeta, rbeta–4.7186 15.6027
ρurho, rrho3.15091
LuLL, rLL–82.9160–88.4363
∇LRscore[0 0][0 7.9145]

For the covariance estimator based on expected Hessian, you need to calculate the expected Hessian analytically. By analytic calculation, the expected value of yi given xi is ρ/βi. Substituting this into Equation 3-2 gives the expected value of the Hessian:

Evaluate this and the resulting covariance estimate numerically as follows:

UEH = [-sum(urho./((ubeta+x).^2)), -sum(1./(ubeta+x)); ...
       -sum(1./(ubeta+x)),-n*UDPsi];
UEstCov3 = -inv(UEH);

Calculate the unrestricted OPG estimator for comparison (lmtest used the restricted OPG estimator):

UG = [-urho./(ubeta+x)+y.*(ubeta+x).^(-2),...
     -log(ubeta+x)-psi(urho)+log(y)];
Uscore = sum(UG)';
UEstCov1 = inv(UG'*UG);

The three estimated covariance matrices, evaluated at the unrestricted maximum, appear in the following table.

OPG (UEstCov1)Observed H–1 (UEstCov2)Expected H–1 (UEstCov3)

These results match those in Greene [28] to a few decimal places. The differences are due to the inaccuracy in finding the maximum likelihood by fminunc with default tolerances, and without providing an analytic gradient in the objective functions nLLF and rnLLF. For higher accuracy, choose smaller TolFun and TolX tolerances, or provide analytic gradients.

Example: Comparing GARCH Models

Test Support for a GARCH(2,1) Model.   Example: Using the Default Model shows that the default GARCH(1,1) model explains most of the variability of the daily returns observations of the Deutschmark/British Pound foreign-exchange rate. This example uses the lratiotest function to determine whether evidence exists to support the use of a GARCH(2,1) model.

The example first fits the Deutschmark/British Pound foreign-exchange rate return series to the default GARCH(1,1) model. It then fits the same series using a GARCH(2,1) model:

The GARCH(1,1) model is the same as the GARCH(2,1) model, except that the parameter G2 is restricted to be 0.

  1. Load the Deutschmark/British Pound foreign-exchange rate data:

    load garchdata
dem2gbp = price2ret(DEM2GBP);

  2. Estimate the GARCH(1,1) model:

    1. Create a GARCH(1,1) default specification structure:

      spec11 = garchset('P',1,'Q',1,'Display','off');

    2. Estimate the model, including the maximized loglikelihood function value, and display the results:

      [coeff11,errors11,LLF11] = garchfit(spec11,dem2gbp);
garchdisp(coeff11,errors11)
Mean: ARMAX(0,0,0); Variance: GARCH(1,1)

Conditional Probability Distribution: Gaussian
Number of Model Parameters Estimated: 4

                               Standard          T     
  Parameter       Value          Error       Statistic 
 -----------   -----------   ------------   -----------
           C    -6.3728e-005   8.379e-005      -0.7606
           K    9.9718e-007    1.2328e-007      8.0890
    GARCH(1)    0.81458        0.015727        51.7956
     ARCH(1)    0.14721        0.013285        11.0812

  3. Estimate the GARCH(2,1) model:

    1. Create a GARCH(2,1) specification structure:

      spec21 = garchset('P',2,'Q',1,'Display','off');

    2. Estimate the GARCH(2,1) model and display the results. Again, calculate the maximized loglikelihood function value.

      [coeff21,errors21,LLF21] = garchfit(spec21,dem2gbp);
garchdisp(coeff21,errors21)
Mean: ARMAX(0,0,0); Variance: GARCH(2,1)

Conditional Probability Distribution: Gaussian
Number of Model Parameters Estimated: 5

                               Standard          T     
  Parameter       Value          Error       Statistic 
 -----------   -----------   ------------   -----------
           C    -4.832e-005    8.4666e-005     -0.5707
           K    1.081e-006     1.4805e-007      7.3014
    GARCH(1)    0.48232        0.11095          4.3472
    GARCH(2)    0.30953        0.10202          3.0341
     ARCH(1)    0.16563        0.016149        10.2564

  4. Perform the Likelihood Ratio Test.

    Since garchfit enforces no boundary constraints during either of the two estimations, you can apply a likelihood ratio test (LR) (see Hamilton [31]). The LR statistic is asymptotically chi-square distributed with degrees of freedom equal to the number of restrictions imposed.

    1. Since the GARCH(1,1) model imposes one restriction, specify one degree of freedom in your call to lratiotest. Test the models at the 0.05 significance level:

      [H,pValue,Stat,CriticalValue] = lratiotest(LLF21,LLF11,...
   1,0.05);
[H,pValue,Stat,CriticalValue]
ans =
    1.0000    0.0210    5.3245    3.8415

      H = 1 indicates that there is enough statistical evidence to support a GARCH(2,1) model.

    2. Alternatively, at the 0.02 significance level:

      [H,pValue,Stat,CriticalValue] = lratiotest(LLF21,LLF11,1,0.02);
[H,pValue,Stat,CriticalValue]
ans =
         0    0.0210    5.3245    5.4119

      H = 0 indicates that there is not enough statistical evidence to support a GARCH(2,1) model.

Testing Multiple Time Series

Multiple Time Series Case Study contains an example that uses vgxvarx to estimate the loglikelihoods of several models, and uses lratiotest to reject some restricted models in favor of an unrestricted model. The calculation appears in Likelihood Ratio Tests.

Back to Top

Akaike and Bayesian Information

You can use Akaike (AIC) and Bayesian (BIC) information criteria to compare alternative models. Information criteria penalize models with additional parameters. Therefore, the AIC and BIC model order selection criteria are based on parsimony (see Box, Jenkins, and Reinsel [10], pages 200-201).

The following example uses the default GARCH(1,1) and GARCH(2,1) models developed in Example: Comparing GARCH Models.

  1. Count the estimated parameters.

    Provide the number of parameters estimated in the model for both AIC and BIC. For the relatively simple models used here, you can just count the number of parameters. The GARCH(2,1) model estimated five parameters:

    • C

    • κ

    • G1

    • G2

    • A1

    The GARCH(1,1) model estimated four parameters:

    • C

    • κ

    • G1

    • A1

      Tip   To count the estimated parameters in more elaborate models, use the function garchcount. garchcount accepts the output specification structure created by garchfit, and returns the number of parameters in the model defined in that structure.

      n21 = garchcount(coeff21)
n21 =
     5
n11 = garchcount(coeff11)
n11 =
     4

  2. Compute the AIC and BIC criteria.

    1. To see the results more precisely, set the numeric format to long:

      format long

        Tip   You can also set the numeric format by selecting File > Preferences > Command Window > Text display from the MATLAB toolbar.

    2. Use the aicbic function to compute the AIC and BIC statistics for the GARCH(2,1) model and the GARCH(1,1) model, and specify the number of observations in the return series:

      [AIC,BIC] = aicbic(LLF21,n21,1974);
[AIC BIC]
ans =
  1.0e+004 *
  -1.596323209285558  -1.593529300675561
[AIC,BIC] = aicbic(LLF11,n11,1974);
[AIC BIC]
ans =
  1.0e+004 *
  -1.595990756407969  -1.593755629519972

      You can use the relative values of the AIC and BIC statistics as guides in the model selection process. In this example, the AIC criterion favors the GARCH(2,1) model, while the BIC criterion favors the GARCH(1,1) default model with fewer parameters. BIC imposes a greater penalty for additional parameters than does AIC. Thus, BIC always provides a given model with a number of parameters no greater than that chosen by AIC.

Back to Top

  

Specification StructuresSpecification Structures Model ConstructionModel Construction

Free Interactive Computational Finance CD

View demos and recorded presentations led by industry experts.

Get the interactive CD

Get trial software

Now On Demand
Network with industry peers and learn the latest applications of the leading software product for computational finance.

 © 1984-2009- The MathWorks, Inc.    -   Site Help   -   Patents   -   Trademarks   -   Privacy Policy   -   Preventing Piracy   -   RSS