nag_correg_robustm_user_varmat (g02hf) calculates an estimate of the asymptotic variance-covariance matrix for the bounded influence regression estimates (M-estimates). It is intended for use with nag_correg_robustm_user (g02hd).

Syntax

[c, wk, ifail] = g02hf(psi, psp, indw, indc, sigma, x, rs, wgt, 'n', n, 'm', m)

[c, wk, ifail] = nag_correg_robustm_user_varmat(psi, psp, indw, indc, sigma, x, rs, wgt, 'n', n, 'm', m)

Description

For a description of bounded influence regression see nag_correg_robustm_user (g02hd). Let

θ

be the regression arguments and let

C

be the asymptotic variance-covariance matrix of

\hat{θ}

. Then for Huber type regression

C = f_{H} {(X^{T} X)}^{- 1} {\hat{σ}}^{2},

where

f_{H} = \frac{1}{n - m} \frac{\sum_{i = 1}^{n} ψ^{2} (r_{i} / \hat{σ})}{{(\frac{1}{n} \sum ψ^{'} (\frac{r_{i}}{\hat{σ}}))}^{2}} κ^{2}

κ^{2} = 1 + \frac{m}{n} \frac{\frac{1}{n} \sum_{i = 1}^{n} {(ψ^{'} (r_{i} / \hat{σ}) - \frac{1}{n} \sum_{i = 1}^{n} ψ^{'} (r_{i} / \hat{σ}))}^{2}}{{(\frac{1}{n} \sum_{i = 1}^{n} ψ^{'} (\frac{r_{i}}{\hat{σ}}))}^{2}},

see Huber (1981) and Marazzi (1987).

For Mallows and Schweppe type regressions,

C

is of the form

{\frac{\hat{σ}}{n}}^{2} S_{1}^{- 1} S_{2} S_{1}^{- 1},

where

S_{1} = \frac{1}{n} X^{T} D X

and

S_{2} = \frac{1}{n} X^{T} P X

D

is a diagonal matrix such that the

i

th element approximates

E (ψ^{'} (r_{i} / (σ w_{i})))

in the Schweppe case and

E (ψ^{'} (r_{i} / σ) w_{i})

in the Mallows case.

P

is a diagonal matrix such that the

i

th element approximates

E (ψ^{2} (r_{i} / (σ w_{i})) w_{i}^{2})

in the Schweppe case and

E (ψ^{2} (r_{i} / σ) w_{i}^{2})

in the Mallows case.

Two approximations are available in nag_correg_robustm_user_varmat (g02hf):

Average over the

r_{i}

\begin{matrix} Schweppe & Mallows \\ D_{i} = (\frac{1}{n} \sum_{j = 1}^{n} ψ^{'} (\frac{r_{j}}{\hat{σ} w_{i}})) w_{i} & D_{i} = (\frac{1}{n} \sum_{j = 1}^{n} ψ^{'} (\frac{r_{j}}{\hat{σ}})) w_{i} \\ P_{i} = (\frac{1}{n} \sum_{j = 1}^{n} ψ^{2} (\frac{r_{j}}{\hat{σ} w_{i}})) w_{i}^{2} & P_{i} = (\frac{1}{n} \sum_{j = 1}^{n} ψ^{2} (\frac{r_{j}}{\hat{σ}})) w_{i}^{2} \end{matrix}

Replace expected value by observed

\begin{matrix} Schweppe & Mallows \\ D_{i} = ψ^{'} (\frac{r_{i}}{\hat{σ} w_{i}}) w_{i} & D_{i} = ψ^{'} (\frac{r_{i}}{\hat{σ}}) w_{i} \\ P_{i} = ψ^{2} (\frac{r_{i}}{\hat{σ} w_{i}}) w_{i}^{2} & P_{i} = ψ^{2} (\frac{r_{i}}{\hat{σ}}) w_{i}^{2} \end{matrix}

See Hampel et al. (1986) and Marazzi (1987).

In all cases

\hat{σ}

is a robust estimate of

σ

nag_correg_robustm_user_varmat (g02hf) is based on routines in ROBETH; see Marazzi (1987).

References

Hampel F R, Ronchetti E M, Rousseeuw P J and Stahel W A (1986) Robust Statistics. The Approach Based on Influence Functions Wiley

Huber P J (1981) Robust Statistics Wiley

Marazzi A (1987) Subroutines for robust and bounded influence regression in ROBETH Cah. Rech. Doc. IUMSP, No. 3 ROB 2 Institut Universitaire de Médecine Sociale et Préventive, Lausanne

Parameters

Compulsory Input Parameters

1: $psi$ – function handle or string containing name of m-file

psi must return the value of the

ψ

function for a given value of its argument.

[result] = psi(t)

Input Parameters

1: $t$ – double scalar: The argument for which psi must be evaluated.

Output Parameters

1: $result$ – double scalar: The value of the function $ψ$ evaluated at t.

2: $psp$ – function handle or string containing name of m-file

psp must return the value of

ψ^{'} (t) = \frac{d}{d t} ψ (t)

for a given value of its argument.

[result] = psp(t)

Input Parameters

1: $t$ – double scalar: The argument for which psp must be evaluated.

Output Parameters

1: $result$ – double scalar: The value of $ψ^{'} (t)$ evaluated at t.

3: $indw$ – int64int32nag_int scalar

The type of regression for which the asymptotic variance-covariance matrix is to be calculated.

$indw < 0$: Mallows type regression.
$indw = 0$: Huber type regression.
$indw > 0$: Schweppe type regression.

4: $indc$ – int64int32nag_int scalar

indw \neq 0

, indc must specify the approximation to be used.

indc = 1

, averaging over residuals.

indc \neq 1

, replacing expected by observed.

indw = 0

, indc is not referenced.

5: $sigma$ – double scalar

The value of

\hat{σ}

, as given by nag_correg_robustm_user (g02hd).

Constraint:

sigma > 0.0

6: $x (ldx, m)$ – double array

ldx, the first dimension of the array, must satisfy the constraint

ldx \geq n

The values of the

X

matrix, i.e., the independent variables.

x (i, j)

must contain the

i j

th element of

X

, for

i = 1, 2, \dots, n

and

j = 1, 2, \dots, m

7: $rs (n)$ – double array

The residuals from the bounded influence regression. These are given by nag_correg_robustm_user (g02hd).

8: $wgt (n)$ – double array

indw \neq 0

, wgt must contain the vector of weights used by the bounded influence regression. These should be used with nag_correg_robustm_user (g02hd).

indw = 0

, wgt is not referenced.

Optional Input Parameters

1: $n$ – int64int32nag_int scalar: Default: the dimension of the arrays rs, wgt and the first dimension of the array x. (An error is raised if these dimensions are not equal.)
$n$ , the number of observations.

Constraint: $n > 1$ .
2: $m$ – int64int32nag_int scalar: Default: the second dimension of the array x.
$m$ , the number of independent variables.

Constraint: $1 \leq m < n$ .

Output Parameters

1: $c (ldc, m)$ – double array: The estimate of the variance-covariance matrix.
2: $wk (m \times (n + m + 1) + 2 \times n)$ – double array: If $indw \neq 0$ , $wk (i)$ , for $i = 1, 2, \dots, n$ , will contain the diagonal elements of the matrix $D$ and $wk (i)$ , for $i = n + 1, \dots, 2 n$ , will contain the diagonal elements of matrix $P$ .
The rest of the array is used as workspace.
3: $ifail$ – int64int32nag_int scalar: $ifail = 0$ unless the function detects an error (see Error Indicators and Warnings).

Error Indicators and Warnings

Errors or warnings detected by the function:

$ifail = 1$

On entry,	$n \leq 1$ ,
or	$m < 1$ ,
or	$n \leq m$ ,
or	$ldc < m$ ,
or	$ldx < n$ .

$ifail = 2$

On entry,

sigma \leq 0.0

$ifail = 3$: If $indw = 0$ then the matrix $X^{T} X$ is either not positive definite, possibly due to rounding errors, or is ill-conditioned.

If $indw \neq 0$ then the matrix $S_{1}$ is singular or almost singular. This may be due to many elements of $D$ being zero.

$ifail = 4$: Either the value of $\frac{1}{n} \sum_{i = 1}^{n} ψ^{'} (\frac{r_{i}}{\hat{σ}}) = 0$ ,

or $κ = 0$ ,

or $\sum_{i = 1}^{n} ψ^{2} (\frac{r_{i}}{\hat{σ}}) = 0$ .

In this situation nag_correg_robustm_user_varmat (g02hf) returns $C$ as ${(X^{T} X)}^{- 1}$ .

$ifail = - 99$: An unexpected error has been triggered by this routine. Please contact NAG.

$ifail = - 399$: Your licence key may have expired or may not have been installed correctly.

$ifail = - 999$: Dynamic memory allocation failed.

Accuracy

In general, the accuracy of the variance-covariance matrix will depend primarily on the accuracy of the results from nag_correg_robustm_user (g02hd).

Further Comments

nag_correg_robustm_user_varmat (g02hf) is only for situations in which

X

has full column rank.

Care has to be taken in the choice of the

ψ

function since if

ψ^{'} (t) = 0

for too wide a range then either the value of

f_{H}

will not exist or too many values of

D_{i}

will be zero and it will not be possible to calculate

C

Example

The asymptotic variance-covariance matrix is calculated for a Schweppe type regression. The values of

X

\hat{σ}

and the residuals and weights are read in. The averaging over residuals approximation is used.

Open in the MATLAB editor: g02hf_example

function g02hf_example


fprintf('g02hf example results\n\n');

x = [1, -1, -1;
     1, -1,  1;
     1,  1, -1;
     1,  1,  1;
     1,  0,  3];

rs  = [0.5643;  -1.1286;   0.5643;  -1.1286;   1.1286];
wgt = [0.4039;   0.5012;   0.4039;   0.5012;   0.3862];

% Estimate variance-covariance matrix
indw  = int64(1);
indc  = int64(1);
sigma = 20.7783;
[c, wk, ifail] = g02hf( ...
                        @psi, @psp, indw, indc, sigma, x, rs, wgt);

disp('Covariance matrix');
disp(c);



function [result] = psi(t)
  if t < -1.5
    result = -1.5;
  elseif abs(t) < 1.5
    result = t;
  else
    result = 1.5;
  end;

function [result] = psp(t)
  if (abs(t) < 1.5)
    result=1;
  else
    result=0;
  end

g02hf example results

Covariance matrix
    0.2070         0   -0.0478
         0    0.2229         0
   -0.0478         0    0.0796

PDF version (NAG web site, 64-bit version, 64-bit version)

Chapter Contents

Chapter Introduction

NAG Toolbox