Forecast-Actuality Significance Test

Introduction

In this notebook we determine the power function for the portmanteau like significance test of Box and Tiao (1976).

We seek to determine the power function for this test for the SIA step intervention with ARIMA disturbances.  We assume then that \(TraditionalForm\`z\_t = ζ + ω\ S\_t\%\((T)\) + \(\(N\_t\)\(.\)\)\)  As before we assume n observations with intervention starting at t = T.  Let m denote the number of post-intervention observations, so m = n - T + 1.

Assuming ω = 0 and taking the forecast origin at T - 1 the lead l forecast may be written,

\(TraditionalForm\`\(z\_\(T - 1\)\)(
        l) = ζ + ψ\_l\ a\_\(T - 1\) + ψ\_\(l + 1\)\ a\_\(T - 2\) +  ... \)

And let (Subscript[e, T - 1])(l) = 
        (Subscript[z, T - 1 + l]) - (Subscript[z, T - 1])(l) be the l step ahead forecast error and let e = ((Subscript[e, T - 1])(
              1),  ... , (Subscript[e, T - 1])(m)).  Then as shown in Box and Tiao (1976, eqn. 4),

\(TraditionalForm\`\(\(Q = \(e\^′\ \(V\^\(-1\)\) 
          e = \(\(e\^′\) π\^′\ π\ e/σ\_a\%2 = \(a\^′\) a/σ\_a\%2\)\)\)\(,\)\)\)

where a = ((Subscript[a, T]),  (Subscript[a, T + 1]),  ... , (Subscript[a, n])) are the innovations, \(TraditionalForm\`\(\(V = \(Var(
            e) = ψ\ \(ψ\^′\) σ\_a\%2\)\)\(,\)\(\ \)\)\)ψ and πare m×m lower triangular matrices.  These matrices may be defined by \(TraditionalForm\`ψ = \((ψ\_\(i - j\))\)\_\(m×m\)\), where the (i, j)-entry is indicated and (Subscript[ψ, k]) is the coefficient of (B^k) in the expansion of (∇^(-d))θ(B)/ϕ(B) = 
        1 + (Subscript[ψ, 1]) B + (Subscript[ψ, 2]) (B^2) +  ... .  Similarly, for π = (Subscript[(((Subscript[π, i - j]))), m×m])if we define (Subscript[π, k]) as the coefficient of (B^k) in the expansion of (∇^d) ϕ(B)/θ(B) = 
        1 + (Subscript[π, 1]) B + (Subscript[π, 2]) (B^2) +  ...

Remark:  Box and Jenkins define (Subscript[π, k]) by (∇^d) ϕ(B)/θ(B) = 
        1 - (Subscript[π, 1]) B - (Subscript[π, 2]) (B^2) -  ...  and this is the definition used in our software.

Under the SIA step model for ω ≠  0,the forecast errors are (Subscript[e, ω]) = ω + e.

(Subscript[Q, ω]) = 
        (Subscript[e, ω]^′) (V^(-1)) (Subscript[e, ω])
=((ω (Subscript[1, m])+e))'(V^(-1))((Subscript[1, m])ω+e)
=((ω (Subscript[1, m])+e))'π (π^′)((Subscript[1, m])ω+e)
=((ω (Subscript[1, m]^′)π/(Subscript[σ, a])+a/(Subscript[σ, a])))'\((\(π\^′\) \(1\_m\) ω/σ\_a + a/σ\_a)\)

Noncentral (χ^2) distribution

In Mathematica the noncentrality parameter, denoted by λ, is defined in Help:

The sum of the squares of n normally distributed random variables with variance (σ^2) = 1 and nonzero means follows a noncentral chi-square distribution NoncentralChiSquareDistribution[n, lambda]. The noncentrality parameter λ is the sum of the squares of the means of the random variables in the sum. Note that in various places in the literature, λ/2 or (Sqrt[λ]) is used as the noncentrality parameter.

According to this definition, \(TraditionalForm\`λ = \((ω\^2/σ\_a\%2)\)\  || \(1\_m\%′\) π\ \( || \^2\)\)

Power Function

Letting δ = ω/σ, the power function can be written,

β(δ) = 
      Pr {X > (Subscript[χ, m]^2)(1 - α) | 
            X~(Subscript[χ, m, λ]^2)}
=1 - 
        F((Subscript[χ, m]^2)(1 - α)),

where F is the noncentral (χ^2)- distribution on m df with noncentrality parameter λ and (Subscript[χ, m]^2)(1 - α) denotes the critical value, ie the 1 - α quantile for the (Subscript[χ, m]^2) distribution.

Computation of the ψ and π matrices

<<Statistics`ContinuousDistributions`

<< d:/math/2003/TSUtil.m

?InvertedCoefficientsARMA

?ImpulseCoefficientsARMA

ψMatrix[m_, phi_, theta_, opts___?Options] := 
      Module[{psi}, 
psi = 
          ImpulseCoefficientsARMA[phi, theta, 
            m - 1]; 
Table[
          PadRight[Reverse[Take[psi, i]], m], {i, m}]]

πMatrix[m_, phi_, theta_, opts___?Options] := 
      Module[{pi}, 
pi = 
          Prepend[-Drop[InvertedCoefficientsARMA[phi, theta, m - 1], 1], 
            1]; 
Table[
          PadRight[Reverse[Take[pi, i]], m], {i, m}]]

ψMatrix[3, {0.5}, {}] // MatrixForm

πMatrix[3, {0.5}, {}] // MatrixForm

Verify that the matrices are inverses.

m = 5; phi = {0.5}; theta = {};

Total[
      Abs[ψMatrix[m, phi, theta] - Inverse[πMatrix[m, phi, theta]]] // 
        Flatten]

Graphical Comparisons

<< d:/math/2004/PowerIA/code.m

?PowerIA

QParameter[m_, phi_, theta_, ω_] := 
      Module[{}, 
(ω^2)*
          Total[Map[Total, 
                Transpose[πMatrix[m, phi, 
                    theta]]]^2]
]

QParameter[20, {0.5}, {}, 1]

PowerQStep[n_, T_, omega_, phi_, theta_] := 
        Module[{m = n - T + 1}, 
lambda = 
            QParameter[m, phi, theta, omega]; 
Xcrit = 
            Quantile[ChiSquareDistribution[m], 0.95]; 
1 - 
            CDF[NoncentralChiSquareDistribution[m, lambda], 
              Xcrit]
];

PowerQStep[120, 100, 1, {0.5}, {}]

PowerQStep[120, 100, -2, {0.5}, {}]

Code Development

phis = {0.0, 0.25, 0.5, 0.75}; j = 3; n = 120; T = 101;

phi = phis[[j]];

htTitle = 0.94; fsZ = 8;

flab = (StyleForm[#1, FontWeight -> "Bold", 
                FontSize -> 
                  fsZ] &) /@ {"δ", "β"};

\(\(pl = 
        StyleForm[
          StringForm[\*""\!\(ϕ\_1\)=``"", 
            NumberForm[phi, {3, 2}]], FontWeight -> "Bold", 
          FontSize -> fsZ];\)\)

Pf1[delta_] := 
      Module[{omega = 
            delta*Sqrt[
                First[TacvfARMA[{phi}, {}, 1]]]}, 
PowerIA[
          n, T, omega, {phi}, {}]
]

Pf2[delta_] := 
      Module[{omega = 
            delta*Sqrt[
                First[TacvfARMA[{phi}, {}, 
                    1]]]}, 
PowerQStep[n, T, 
          omega, {phi}, {}]
]

Plot[{Pf1[delta], Pf2[delta]}, {delta, 0, 
          2}, 
Prolog -> {Text[
              pl, {0, htTitle}, {0, 
                0}]}, 
AxesOrigin -> {0, 0}, 
        Frame -> True, FrameLabel -> flab, 
        RotateLabel -> False, 
        FrameTicks -> {{0, 2}, {0.25, 0.75}, None, None}, 
        Axes ->  False, TextStyle -> {FontSize -> fsZ}, 
        PlotStyle -> {Dashing[{}], Dashing[{0.03, 0.02}]}];

Final Code

phis = {0.0, 0.4, 0.8}; n = 120; T = 101;

htTitle = 0.94; fsZ = 8;

flab = (StyleForm[#1, FontWeight -> "Bold", 
                FontSize -> 
                  fsZ] &) /@ {"δ", "β"};

Pf1[delta_] := 
      Module[{}, 
omega = 
          delta*Sqrt[
              First[TacvfARMA[{phi}, {}, 1]]]; 
PowerIA[n, 
          T, omega, {phi}, {}]
]

Pf2[delta_] := 
      Module[{}, 
omega = 
          delta*Sqrt[
              First[TacvfARMA[{phi}, {}, 1]]]; 
PowerQStep[
          n, T, omega, {phi}, {}]
]

g = Array[0 &, Length[phis]];
\(Do[
phi = 
        phis[[
          j]]; 
pl = 
        StyleForm[
          StringForm[\*""\!\(ϕ\_1\)=``"", 
            NumberForm[phi, {3, 2}]], FontWeight -> "Bold", 
          FontSize -> fsZ]; 
g[[
          j]] = 
        Plot[{Pf1[delta], Pf2[delta]}, {delta, 0, 
            2}, 
Prolog -> {Text[
                pl, {1.5, 0.15}, {0, 
                  0}]}, 
AxesOrigin -> {0, 0}, 
          Frame -> True, FrameLabel -> flab, 
          RotateLabel -> False, 
          FrameTicks -> {{0, 1, 2}, {0.25, 0.75}, None, None}, 
          Axes -> \ False, TextStyle -> {FontSize -> fsZ}, 
          PlotStyle -> {Dashing[{}], 
              Dashing[{0.03, 0.02}]}], 
{j, 
        Length[phis]}]\)

Graphics Array

gOUT = 
        Show[GraphicsArray[g, GraphicsSpacing -> -0.0, 
            FrameTicks -> {{0.25, 0.5, 0.75}, {-2, 2}, None, 
                None}]];

Export Graph

Export["d:/math/2004/PowerIA/Simulation/Q.eps", 
        gOUT, "EPS"];

Difference Graph

g2 = Array[0 &, Length[phis]];
\(Do[
phi = 
        phis[[
          j]]; 
pl = 
        StyleForm[
          StringForm[\*""\!\(ϕ\_1\)=``"", 
            NumberForm[phi, {3, 2}]], FontWeight -> "Bold", 
          FontSize -> fsZ]; 
g2[[
          j]] = 
        Plot[{Pf1[delta] - Pf2[delta]}, {delta, 0, 
            2}, 
Prolog -> {Text[
                pl, {1.5, 0.15}, {0, 
                  0}]}, 
AxesOrigin -> {0, 0}, 
          Frame -> True, FrameLabel -> flab, 
          RotateLabel -> False, 
          FrameTicks -> {{0, 1, 2}, {0.25, 0.75}, None, None}, 
          Axes -> \ False, TextStyle -> {FontSize -> fsZ}, 
          PlotStyle -> {Dashing[{}], 
              Dashing[{0.03, 0.02}]}], 
{j, 
        Length[phis]}]\)

Graphics Array

gOUT = 
        Show[GraphicsArray[g2, GraphicsSpacing -> -0.0, 
            FrameTicks -> {{0.25, 0.5, 0.75}, {-2, 2}, None, 
                None}]];

Tabulate Power Functions for Input to S-Plus Trellis

i =. ; j =. ; deltas = Range[0.0, 2, 0.01]; 
    g = Array[0 &, Length[phis]];
Do[
phi = 
        phis[[
          j]]; 
g[[
          j]] = 
        Table[{Pf1[deltas[[i]]], 
            Pf2[deltas[[i]]]}, {i, 
            Length[deltas]}], 
{j, 
        Length[phis]}]

g = Flatten[g, 1];

Reference

Comparison of Forecast and Actuality G. E. P. Box, G. C. Tiao Applied Statistics, Vow. 25, No. 3. (1976), pop. 195-200.