RobCov

robCov computes covariance matrix of robust regression coefficients

Syntax

  • out=RobCov(X,scaledres,scaleest)example
  • out=RobCov(X,scaledres,scaleest,Name,Value)example

Description

Under some regularity conditions, robust (S and MM) estimates are asymptotically normal, thereby allowing for Wald-type tests and confidence intervals. The covariance matrix of the estimated parameters \[ cov(\hat \beta)= q^2 \times \sigma^2 \times v \times V_X^{-1} \]

consists of four parts:

1) $q$ a correction factor for the scale estimate;

2) $\sigma$ the scale parameter.

3) $v$ a correction factor depending on the $\psi$ function which is used;

4) $V_X$= a matrix part. For OLS $V_X=X'X$. Given that in robust regression we give a weight to each observation, the matrix $X'X$ should be replaced by something like $X'WX$, where $W$ is a diagonal matrix containing the weights assigned to each observation.

The purpose of this function is to provide the user with different options for the estimate of $cov(\hat \beta)$ where $\hat \beta$ is a vector of regression coefficients obtained using S or MM estimation and a particular $\rho$ function.

example

out =RobCov(X, scaledres, scaleest) Compare the 5 estimates of cov matrix.

example

out =RobCov(X, scaledres, scaleest, Name, Value) Compare t stat from S and MM estimator.

Examples

expand all

  • Compare the 5 estimates of cov matrix.
  • n=200;
    p=3;
    state1=123456;
    randn('state', state1);
    X=randn(n,p);
    y=randn(n,1);
    kk=10;
    ycont = y;
    ycont(1:kk)=ycont(1:kk)+7;
    [outS]=Sreg(ycont,X);
    rhofunc='optimal';
    bdp=0.5;
    out=RobCov(X,outS.residuals,outS.scale);
    disp('Compare 5 estimates of cov(\hat beta)')
    disp(out.covrob)
    disp('--------')
    disp(out.covrob1)
    disp('--------')
    disp(out.covrob2)
    disp('--------')
    disp(out.covrob3)
    disp('--------')
    disp(out.covrob4)

  • Compare t stat from S and MM estimator.
  • rhofunc='optimal';
    bdp=0.5;
    out=RobCov(X,outS.residuals,outS.scale,'rhofunc',rhofunc,'bdp',0.5);
    covrobS=out.covrob;
    covrobS1=out.covrob1;
    covrobS2=out.covrob2;
    covrobS3=out.covrob3;
    covrobS4=out.covrob4;
    % Compute robust S t-statistics
    tstatS=outS.beta./sqrt(diag(covrobS));
    tstatS1=outS.beta./sqrt(diag(covrobS1));
    tstatS2=outS.beta./sqrt(diag(covrobS2));
    tstatS3=outS.beta./sqrt(diag(covrobS3));
    tstatS4=outS.beta./sqrt(diag(covrobS4));
    eff=0.95;
    outMM=MMregcore(ycont,X,outS.beta,outS.scale);
    out=RobCov(X,outMM.residuals,outS.scale,'rhofunc',rhofunc,'eff',eff);
    covrobMM=out.covrob;
    covrobMM1=out.covrob1;
    covrobMM2=out.covrob2;
    covrobMM3=out.covrob3;
    covrobMM4=out.covrob4;
    tstatMM=outMM.beta./sqrt(diag(covrobMM));
    tstatMM1=outMM.beta./sqrt(diag(covrobMM1));
    tstatMM2=outMM.beta./sqrt(diag(covrobMM2));
    tstatMM3=outMM.beta./sqrt(diag(covrobMM3));
    tstatMM4=outMM.beta./sqrt(diag(covrobMM4));
    disp('tstat from S')
    disp([tstatS tstatS1 tstatS2 tstatS3 tstatS4])
    disp('--------')
    disp('tstat from MM')
    disp([tstatMM tstatMM1 tstatMM2 tstatMM3 tstatMM4])
    qhat=out.q;
    disp('tstat from MM after correction for sigma')
    disp([tstatMM/qhat tstatMM1/qhat tstatMM2/qhat tstatMM3/qhat tstatMM4/qhat])

    Related Examples

    expand all

  • An example of the need of using covrobc.
  • XX is a 50-by-2 matrix generated by standard random normal.

    % Estimate beta using Sreg and compare the four estimators of 
    % cov (\hat \beta).
    XX=[0.8115    1.3104
    -0.4477    0.8110
    0.9438    1.9121
    -0.2346    0.7058
    1.8364    3.1318
    0.2911    0.7837
    0.3669    1.7571
    0.4560    1.0876
    -1.0831    0.0679
    0.7018    1.0923
    1.6976    2.7744
    -0.0624    1.2015
    0.2016    0.7747
    -0.1542    1.3095
    -0.9416   -0.5251
    0.4873    0.8116
    -0.3561    0.9261
    0.5185    2.1219
    -0.7674   -0.0514
    0.3360    1.8620
    -0.2083    1.3796
    -0.1700    0.9053
    1.9152    2.4369
    -0.9912   -0.6885
    0.5939    2.0449
    -1.1842   -0.8506
    1.1375    1.6184
    -0.2482    1.2104
    0.7936    2.2313
    -0.2744    1.5259
    1.1455    2.2739
    -0.0637    1.2618
    -0.5054    0.0536
    0.7282    1.9411
    -0.0703    0.8090
    -1.8186   -0.8725
    -2.0621   -0.5299
    -1.7750   -0.7765
    -0.8553   -0.3311
    1.8966    2.9885
    -0.1076    0.5106
    0.9055    1.1283
    0.3086    1.5491
    -0.0487    0.4964
    -1.6901   -0.7348
    -0.7253    0.9436
    0.2124    0.6549
    -0.2748   -0.2009
    -1.2804   -0.7036
    0.3627    1.3027];
    % Define X and y;
    y=XX(:,2);
    X=XX(:,1);
    outS=Sreg(y,X,'rhofunc','bisquare','bdp',0.5);
    outRC=RobCov(X,outS.residuals,1,'rhofunc','bisquare','bdp',0.5);
    disp('covrob')
    disp(outRC.covrob)
    disp('covrob1')
    disp(outRC.covrob1)
    disp('covrob2')
    disp(outRC.covrob2)
    disp('covrob3')
    disp(outRC.covrob3)
    disp('covrob4')
    disp(outRC.covrob4)
    disp('covrobc')
    disp(outRC.covrobc)
    % plot X and y and add LS and robust regression line
    plot(X,y,'o')
    xlabel('X')
    ylabel('y')
    hLS=lsline;
    hLS.Color='r';
    hLS.LineStyle = ':';
    hLS.LineWidth = 2;
    hROB=refline([outS.beta(2) outS.beta(1)]);
    hROB.LineWidth = 2;
    legend({'' 'LS','Robust'},'Location','Best')
    Total estimated time to complete S estimate:  1.74 seconds 
    covrob
        0.6428   -0.0036
       -0.0036    0.7251
    
    covrob1
       29.6408   -0.4576
       -0.4576   29.4429
    
    covrob2
       29.6345   -0.1655
       -0.1655   33.4256
    
    covrob3
        0.2392   -0.0037
       -0.0037    0.2376
    
    covrob4
       1.0e-03 *
    
        0.2844   -0.0069
       -0.0069    0.2488
    
    covrobc
        0.6428   -0.0036
       -0.0036    0.7251
    
    
    Click here for the graphical output of this example (link to Ro.S.A. website)

    Input Arguments

    expand all

    X — Data matrix of explanatory variables (also called 'regressors') of dimension (n x p-1). Rows of X represent observations, and columns represent variables.

    Data Types: single | double

    scaledres — Scaled residuals.Vector. n-times-1 vector containing scaled residuals $r_i/\hat \sigma$.

    Data Types: single | double

    scaleest — robust estimate of the scale. Scalar.

    Robust estimate of sigma ($\hat \sigma$).

    Data Types: single | double

    Name-Value Pair Arguments

    Specify optional comma-separated pairs of Name,Value arguments. Name is the argument name and Value is the corresponding value. Name must appear inside single quotes (' '). You can specify several name and value pair arguments in any order as Name1,Value1,...,NameN,ValueN.

    Example: 'intercept',1 , 'eff',0.99 , 'intercept',1 , 'bdp',0.4 , 'rhofunc','optimal' , 'rhofuncparam',5

    intercept —Indicator for constant term.scalar.

    If 1, a model with constant term will be fitted (default), else no constant term will be included.

    Example: 'intercept',1

    Data Types: double

    eff —nominal efficiency.scalar.

    Scalar defining nominal efficiency (i.e. a number between 0.5 and 0.99). The default value is 0.95 Asymptotic nominal efficiency is:

    $(\int \psi' d\Phi)^2 / (\psi^2 d\Phi)$

    Example: 'eff',0.99

    Data Types: double

    intercept —Indicator for constant term.scalar.

    If 1, a model with constant term will be fitted (default), else no constant term will be included.

    Example: 'intercept',1

    Data Types: double

    bdp —breakdown point.scalar.

    It measures the fraction of outliers the algorithm should resist. In this case any value greater than 0 but smaller or equal than 0.5 will do fine.

    Note that given bdp nominal efficiency is automatically determined.

    REMARK: just one between bdp and eff must be specified. If both of them are specified an error is produced. If both of them are not specified the defulat is tu use the tuning constant associated to a nominal efficiency of 0.95.

    Example: 'bdp',0.4

    Data Types: double

    rhofunc —rho function.string.

    String which specifies the rho function which must be used to weight the residuals. Possible values are 'bisquare' 'optimal' 'hyperbolic' 'hampel'.

    'bisquare' uses Tukey's $\rho$ and $\psi$ functions.

    See TBrho and TBpsi.

    'optimal' uses optimal $\rho$ and $\psi$ functions.

    See OPTrho and OPTpsi.

    'hyperbolic' uses hyperbolic $\rho$ and $\psi$ functions.

    See HYPrho and HYPpsi.

    'hampel' uses Hampel $\rho$ and $\psi$ functions.

    See HArho and HApsi.

    The default is bisquare

    Example: 'rhofunc','optimal'

    Data Types: double

    rhofuncparam —Additional parameters for the specified rho function.scalar | vector.

    For hyperbolic rho function it is possible to set up the value of k = sup CVC (the default value of k is 4.5).

    For Hampel rho function it is possible to define parameters a, b and c (the default values are a=2, b=4, c=8)

    Example: 'rhofuncparam',5

    Data Types: single | double

    Output Arguments

    expand all

    out — description Structure

    A structure containing the following fields

    Value Description
    covrob

    p-times-p (if intercept is 1 else is (p-1)-by-(p-1)) matrix containing asymptotic variance covariance matrix of regression coefficients. covrob implements equation (4.49) of p. 101 of Maronna et al. (2006) namely:

    \[ \mbox{covrob} = cov( \hat \beta) = \hat \sigma^2 \hat v (X'X)^{-1} \] where \[ \hat v = \frac{n}{n-p} n\frac{\sum_{i=1}^n \psi(r_i/\hat \sigma)^2}{\sum_{i=1}^n \psi'(r_i/\hat \sigma)^2} \]

    covrob1

    p-times-p (if intercept is 1 else is (p-1)-by-(p-1)) matrix containing asymptotic variance covariance matrix of regression coefficients. covrob1 implements equation (7.81) of p. 171 of Huber and Ronchetti (2009) with $(X'X)^{-1}$ replaced by $(X' W X)^{-1}$ namely:

    \[ \mbox{covrob1} = K^2 \hat v (X' W X)^{-1}; \] where $K=1+p n \frac{var(\psi' (r/\hat \sigma))}{ \left[ (\sum_{i=1}^n \psi'(r_i/\hat \sigma)) \right]^2}$ ;

    covrob2

    p-times-p (if intercept is 1 else is (p-1)-by-(p-1)) matrix containing asymptotic variance covariance matrix of regression coefficients. covrob1 implements equation (7.81) of p. 171 of Huber and Ronchetti (2009) with $X'X$ and $K^2$ namely:

    \[ \mbox{covrob2} = K^2 \hat v (X' X)^{-1}; \] where $K=1+p n \frac{var(\psi' (r/\hat \sigma))}{ (\sum_{i=1}^n \psi'(r_i/\hat \sigma))^2}$ ;

    covrob3

    p-times-p (if intercept is 1 else is (p-1)-by-(p-1)) matrix containing asymptotic variance covariance matrix of regression coefficients. covrob implements equation (7.82) of p. 171 of of Huber and Ronchetti (2009).

    namely:

    \[ \mbox{covrob3} = K \hat v (X' W X)^{-1}; \] where $K=1+p n \frac{var(\psi' (r/\hat \sigma))}{ (\sum_{i=1}^n \psi'(r_i/\hat \sigma))^2}$ ;

    covrob4

    p-times-p (if intercept is 1 else is (p-1)-by-(p-1)) matrix containing asymptotic variance covariance matrix of regression coefficients. covrob implements equation (7.83) of p. 171 of of Huber and Ronchetti (2009).

    namely:

    \[ \mbox{covrob4} = \frac{1}{n-p} K^{-1} \sum_{i=1}^n (\psi(r_i/\hat \sigma))^2 (X' W X)^{-1} X'X (X' W X)^{-1}; \]

    covrobc

    p-times-p (if intercept is 1 else is (p-1)-by-(p-1)) matrix containing best variance covariance matrix of regression coefficients. covrobc checks the coefficient of variation of the derivative of function $\psi(r_i/\hat \sigma)$. If it is greater than the threshold proposed in Salini, et al. (2019), it uses covrob2 else covrob.

    q

    scalar. Correction for scale estimate (see Maronna and Yohai CSDA 2010). It is defined as \[ q=1+\frac{p}{2n} \frac{a}{b \times c} \] where \[ a = \frac{1}{n} \sum_{i=1}^n (\psi(r_i/\hat \sigma))^2 \] \[ b = \frac{1}{n} \sum_{i=1}^n (\psi'(r_i/\hat \sigma))^2 \] \[ c = \frac{1}{n} \sum_{i=1}^n (\psi(r_i/\hat \sigma)) r_i/\hat \sigma) \]

    References

    Maronna, R.A., Martin D. and Yohai V.J. (2006), "Robust Statistics, Theory and Methods", Wiley, New York.

    Huber, P.J. and Ronchetti, E.M. (2009), "Robust Statistics, 2nd Edition", Wiley.

    Maronna, R.A., and Yohai V.J. (2010), Correcting MM estimates for fat data sets, "Computational Statistics and Data Analysis", Vol. 54, pp. 3168-3173.

    Koller, M. and W. A. Stahel (2011), Sharpening wald-type inference in robust regression for small samples, "Computational Statistics & Data Analysis", Vol. 55, pp. 2504-2515.

    Croux, C., Dhaene G., and Hoorelbeke D. (2003), Robust standard errors for robust estimators. Technical report, Dept. of Applied Economics, KU Leuven.

    Salini, S., Laurini, F., Morelli, G., Riani M. and Cerioli A. (2019), Covariance matrices of robust estimators in regression, submitted.

    See Also

    | |

    This page has been automatically generated by our routine publishFS