Bias of Burg Estimate in AR(1)

Introduction and Purpose

In this notebook the asymptotic mean and variance for the Burg estimate \(TraditionalForm\`ϕ\&^\) in an AR(1) is derived.  It is shown that \(TraditionalForm\`E(ϕ\&^) = \(-2\)\ ϕ/n + \(O\_p\)(
            n\^\(-2\))\).  As a check on our computations the asymptotic variance which is well-known to be given by \(TraditionalForm\`Var(ϕ\&^) = \((1 - ϕ\^2)\)/n\) is also dervied.

The bias for several estimators, unknown mean case, of the AR(1) are compared in the table below.

({{Yule - Walker, -(1 + 4  ϕ)/n}, {Least - Squares, -(1 + 3) ϕ/n}})

Initializations

Acvf AR(1)

Let γ[i,j] denote the covariance of (Subscript[x, i]) and (Subscript[x, j]).  So for the AR(1), γ[i, j] = (ρ^(|i - j|)).

γ[i_, 
          j_] = (ρ^Abs)[i - j]/(1 - (ρ^2));

Isserlis Formula

Isserlis (1918) has derived a formula for the product-moment coefficient of any order of a normal distribution frequency in any number of variables.

Isserlis[{i_, j_}] := γ[i, j];

Isserlis[a_] := 
      Module[{k, na = Length[a]}, 
        Sum[Isserlis[{a[[1]], 
                a[[1 + k]]}] 
            Isserlis[Drop1k[a, k]], {k, na - 1}]]

Drop1k[a_, k_] := Drop[Drop[a, {1}], {k}];

Let "a__" reprerent an index vector {(Subscript[i, 1]),  (Subscript[i, 2]),  ... ,  (Subscript[i, k])} of a time series product (Subscript[x, (Subscript[i, 1])]) 
            (Subscript[x, (Subscript[i, 2])]) ...  (Subscript[x, (Subscript[i, k])]).    Isserlis[a_] will compute the expectation of that time series prodect, E ((Subscript[x, (Subscript[i, 1])]) (Subscript[x, (Subscript[i, 2])]) ...  
            (Subscript[x, (Subscript[i, k])])).  For example, Isserlis[{1, 2, 3, 4}] = 
        E ((Subscript[x, 1]) (Subscript[x, 2]) (Subscript[x, 3]) (Subscript[x, 4])).

Isserlis[{t, t + 1, t + 2, t + 3}] // Simplify

Isserlis[{1, 2, 3, 4}] // Simplify

CoeQ.  Test if it is a coefficient and not a RV

CoeQ.  Coefficients of (ρ^k)

CoeQ[c_Integer (ρ^k)_] := True;

CoeQ[(ρ^k)_] := True;

CoeQ[c_Integer] := True;

CoeQ[c_Integer ρ] := True;

CoeQ[ρ] := True;

IsserlisProduct

For convenience, Isserlis[a__] is modified to be IsserlisProduct[a__] that allows to directly take the product of time series (Subscript[x, (Subscript[i, 1])]) (Subscript[x, (Subscript[i, 2])]) ...  
        (Subscript[x, (Subscript[i, k])]) as the input.

Pow[ω_ a__]:=ω Pow[Times[a]]/;CoeQ[ω];

HoldPattern[Pow[Times[a__]]]:=Flatten[Map[Pow,List[a]]];

Pow[x[k_]^p_]:=Table[k,{p}];

Pow[x[k_]]:={k};

Protect[Pow];

IsserlisProduct[ω_ a__] := ω IsserlisProduct[Times[a]] /; 
        CoeQ[ω]

IsserlisProduct[a__] := Isserlis[Pow[a]]

IsserlisProduct[x[t] x[t + 1] x[t + 2] x[t + 3]] // Simplify

Burg Estimate

The Burg AR(1) estimate may be written (Overscript[ρ, ∧]) = A/B, where

A=(Overscript[((∑_(t = 2))), n])2  x[t] x[t + 1]

and

B=(Overscript[((∑_(t = 2))), n])x[t]^2+x[t + 1]^2

So we define the terms,

Bt[t_] := x[t]^2 + x[t + 1]^2;

At[t_] := 2  x[t] x[t + 1];

Bias Computation

E ((Overscript[ρ,  ∧]) ) can be expanded  to order (n^(-1)) as

E ((Overscript[ρ, ∧]))=((Subscript[μ, A])/(Subscript[μ, B]))+\(\(1\/\((μ\_B)\)\^2\) {E\ \((B\^2)\)\  - \((μ\_B)\)\^2}\ μ\_A\/μ\_B\)-((E(A B) - (Subscript[μ, A]) (Subscript[μ, B]))/((((Subscript[μ, B])))^2))+O((1/(n^2))).

Normal[Series[A/B - mA/mB, {A, mA, 1}, {B, mB, 1}]]

Normal[Series[A/B - mA/mB, {A, mA, 1}, {B, mB, 1}]] // 
      FullSimplify

To compute the expectation of B, (Subscript[μ, B]):

μB = (∑_(t = 1)^(n - 1))Map[IsserlisProduct,  
          Expand[Bt[t]]]

To compute the expectation of A, (Subscript[μ, A]):

μA = (∑_(t = 1)^(n - 1))IsserlisProduct[At[t]]

To compute the expectation of (B^2), E ((B^2)):

μBsq = (∑_(t = 1)^(n - 1))Map[IsserlisProduct,  
              Expand[Bt[t]^2]] + 
          2 (∑_(t = 1)^(n - 2))(∑_(j = 1)^(n - t))Map[
                  IsserlisProduct, Expand[Bt[t] Bt[t + j]]] // 
        Simplify

To compute the expectation of (A^2), E ((A^2)):

μAsq = (∑_(t = 1)^(n - 1))IsserlisProduct[At[t]^2] + 
          2 (∑_(t = 1)^(n - 2))(∑_(j = 1)^(n - t))IsserlisProduct[At[t] At[t + j]] // Simplify

To compute the expectation of A B, E (A B):

μAB = (∑_(t = 1)^(n - 1))Map[IsserlisProduct, 
              Expand[Bt[t] At[
                    t]]] +  (∑_(t = 1)^(n - 2))(∑_(j = 1)^(n - t)) 
              Map[IsserlisProduct, 
                Expand[Bt[t] 
                    At[t + j]]] +  (∑_(t = 1)^(n - 2))(∑_(j = 1)^(n - t)) 
              Map[IsserlisProduct, Expand[Bt[t + j] At[t]]] // 
        Simplify

Then ((Subscript[μ, Y])/(Subscript[μ, X])) is

Mratio = μ(A/μ)B

Therefore, the bias of \(TraditionalForm\`\(ϕ\&^\)\) to order (n^(-1)) can be derived as

E ((Overscript[ρ, ∧]) - ρ)≐\(\(\(\(\ \)\(μ\_A\)\)\/\((μ\_B)\)\^3\) {E\ \((B\^2)\)\  - \((μ\_B)\)\^2}\)-((E(A B) - (Subscript[μ, A]) (Subscript[μ, B]))/((((Subscript[μ, B])))^2))

BurgBias = 
Expand[μ(A/(((μB))^3)) (μBsq - (((μB))^2)) - (μAB - (μA) (μB))/(((μB))^2) // Simplify]

The bias of Burg estimate for AR(1) to order (n^(-1)) is

BurgBias = Apart[BurgBias, n]

Limit[n*BurgBias, n -> ∞, 
      Assumptions -> -1 < ρ < 0]

Bias Computation, Trying another expression

Normal[Series[A/B - mA/mB, {A, mA, 1}, {B, mB, 1}]] // 
      FullSimplify

BurgBiasSimpler = 
Expand[μ(A/(((μB))^3)) (μBsq - (((μB))^2)) - (μAB - (μA) (μB))/(((μB))^2) // Simplify]

BurgBiasSimpler = 
Expand[-(μAB - 
                  2 μA μB  + μA μB)/(((μB))^2) // 
          Simplify]

BurgBiasSimpler = Apart[BurgBias, n]

Limit[n*BurgBiasSimpler, n -> ∞, 
      Assumptions -> -1 < ρ < 0]

It appears that only μAB and μB are needed.  It appears the other terms vanish.

BurgBiasEvenSimpler = 
Expand[-(μAB)/(((μB))^2) // Simplify]

BurgBiasSimpler = Apart[BurgBias, n]

Limit[n*BurgBiasEvenSimpler, n -> ∞, 
      Assumptions -> -1 < ρ < 0]

Variance Computation

The variance of (Overscript[ρ,  ∧])  to order (n^(-1)) can be derived as

\(TraditionalForm\`Var(\(ρ\& ∧ \)) = \(\(\((E\ \((X\^2)\) - \((μ\_X)\)\^2)\) \((μ\_\(\(\ \)\(Y\)\))\)\^2\/\((μ\_X)\)\^4\)\(+\)\(\((E\ \((Y\^2)\) - \((μ\_\(\(\ \)\(Y\)\))\)\^2)\)\/\((μ\_X)\)\^2\)\(\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \)\)\)

\(TraditionalForm\`\(\(\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \)\(\(\(-2\) \((E(X\ Y) - μ\_X\ μ\_Y)\) μ\_Y\/\((μ\_X)\)\^3 + 
      O(n\^\(-2\))\)\(\ \)\(,\)\(\ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \ \)\)\)\)

BurgVariance = 
      Expand[(uxx - (((ux))^2)) ((((uy))^2)/(((ux))^4)) + (((uyy - (((uy))^2)))/(((ux))^2)) - 
            2 (uxy - (ux) (uy)) (uy/(((ux))^3)) // 
          Simplify]

The variance of Burg estimate for AR(1) to order (n^(-1)) is

BV = Apart[BurgVariance, n]

Limit[n*BV, n -> ∞, 
      Assumptions -> -1 < ρ < 0]