This paper gives an alternative derivation of a Monte Carlo method that has been used to study robust estimators. Extensions of the technique to the regression case are also considered and some computational points are briefly mentioned.

Paul W. Holland

English

Research Papers in Economics

NBER WORKINGPAPER SERIES
MONTECARLO FOR ROBUST REGRESSION:
THE SWINDLE UNMASKED
PaulW..Holland*
WorkingPaperNo. 10
COMPUTER RESEARCHCENTER FOR ECONOMICSAND MANAGEMENTSCIENCE National Bureau ofEconomicResearch,Inc.
575Technology Square
Cambridge, Massachusetts02139
September 1973
Preliminary: not for quotation
NBER working papers are distributedinformally and in limited numbers
for comments only. They should not bequoted without written permission.
Thisreport has not undergone the review accorded officialNBER publications; in particular, it has not yet beensubmitted for approval by the Board of Directors.
*NBER Computer Research Center. Researchsupported inpart by National ScienceFoundation Grant GJ-1154X2 to the NationalBureau cf Economic Research, Inc.Abstract
Givesan alternative derivation of a tlonte Carlo method that has been
used to study robust estimators. Extensions of thetechnique to the regressthn
case are also considered and somecomputational points are briefly mentioned.
(
(1.Introduction
In this paper, wediscussa method for achieving moreaccuracy from a
Monte Carlo study than is possible from
simple random sampling.. Such Monte
Carlo "swindles" are important in thelarge scale use of Monte Carlo studies.
The particular method we discuss herehas been described before in Relies[1970]
and Andrews etal.[1972],but their approaches are somewhat different from
the one we employ. A deeperunderstanding of the method and its properties
is gained by having alternative derivationsavailable.
The particular problem weconsideris the following. Webeginwith the
familiar linear regression problem
y =X8+e
(1-1)
where y is Nxl, X is Nxp,8is pxl and e is Nxl. We furthermoreassume that
the components of e, e., are independentand identically distributed random
variables with common density
:•f() (1-2)
where f is assumed to be symmetric abou0, i.e., f(-.x) =f(x).The linear
regression problem is to estimate 8 usingy and X. Let 8 denote a generic
estimator of 8. Sometimes its dependenceon y will be denoted by (y).
There are two notions of invariance that willbe important in the rest
of this paper.
Scale Invariance
An estimator 8(y) is said to be scale invariantif
8(c y) C(y)
(1-3)
for any constant, c.-2—
Rerjrssion Invariance
An estiiator (y) issaidto be regression (on X) invariant if
(y+.x1)(y) + y (1-4)
foranypxl vector, •y.
We shall restrict our discussion to estimators, ,whichare both regres-
sion and scale invariant. The problem of main concern is to study
•Cov) = - - (1-5)
Howeverbecause we have restricted attention to.regression and scale invariant
estimators we have
E(y) - - ))T2E((Y - X)T) =
• (1-6)
Thus we.may assumewithout loss of generality that 0 and 2
In order to compute Cov0(B), we must often resort to a Monte Carlo study.
The swindle we will consider is designed for such an investigation. When
p1 and x.1, the regression problem reduces to the "location" problem in
which we are estimating the center of a symmetric distribution. In
thelocation case (1-5) reduces to the variance of .
We shall divide our discussion of the swindle into fourcases: Location
withGaussian errors, Location with Gauss/independent errors, Regression with
Gaussian errors and Regression with Gauss/Independent errors.2. Location with Gaussian Errors
In this case we are concerned with computing
Var0() =E0(2) (2-1).
for a general location and scale invariant estimator underthe assumption that
f is the unit Gaussian density. However, when f is the Gaussiardistribution,
we know what the best location (regression) and scale invariantestimator is ——
ourfriend ,thesample average. Hence, instead of computing the variance of
,we shalltry to compute the excess of the variance of over that of ,.We
now derive the important formula that allows us todo this.
Theorem 1: Ifis ylocationand scale invariant estimator, then under
unit Gaussian errors we have
Var(s) = Var() + E( -)2 (2-2)
p.
= Var()+E((5: (2-3)
where 2 is the usualunbiased estimate of o.
Proof:We begin with (2—2) and then derive (2-3) from it. We have
((y))2=((y-) + )2 = ((y-))2 + 2 + ()2
(2-4)
Hence
Var(s)E0(2)E0()2 +2E0( (y -)) +
E0(B(y
(2—5)
But andy -areindependent so thatand (y -) areindependent and
hence .
E0((y -)) = 0-4-
Also E0G)2Var() and (y - =-so that Vr()=Var() +E(-
which was to be proved. To prove (2—3) we need to show that
E( - = (2-6)
This is done as follows
• E( -•= E(S(Y))2) (2-7)
FtS is independent of so that S2 is independent of(Y)') and
hence(2—7)equals
E(s2) E( Y)2)2 E(s
)2
(2-8)
Sincewehave assumed G2=1,(2-6) follows immediately. OED
It should be noted that (2-2) onlyrequiresto be location invariant
while (2-3) requires both location and scale invariance.
it is quite simple to use Theorem 1 to get a Monte Carlo swindle for the
variance of a location and scale invariant estimator, B. If weweregoing to
use naiverandomsampling to estimate Var(s) wewoulddraw repeated independent
samples of size N, y =(y1,...,YN)T,compute(y) for each sample and
average the value of 2 over the replications. But from (2—2) weseethat we
may also estimate Var(s) by computing (y) and .fromeach sample, averaging
the value of (— )2over the replications and then adding N1 =Var()to
the result, i.e.,
estimated Var(s) =N +Monte Carlo Average (—
(2-9)
In order to remove the first order effect of N when looking at different
sample sizes it is customary to estimate N Var(B) instead This leads to— 5-.
estimated(Var(s)) =1+rI(MonteCarlo Average (- )2)
(2-10)
What is the nature of this swindle? The mainpoint is that it concentrates
the slow Monte Carlo convergence on theexcess of the Var(B) over Var() rather
than allowing this to effect all of Var(S). This implies that if B is a
very efficient estimator (has a small excess of variance) then the swindle will
be more effective than if 5 has a largeexcess. This is because the smaller
the excess, the smaller the percentage of theVar(S) that is estimated by
Monte Carlo averaging. From (2—4) and (2—5)weseethat the swindle utilizes
the theory of the Gaussian distribution toget exact resu1s for twoofthree
pieces of E(2); theseare
E( (y -)) = 0
(2-11)
andE(5)2 =f11.
.
(2-12)
Monte Carlo is then used to estimate the thirdpiece.
We nay use (2-3) to get more of a swindle via theestimate:
estimated N(Var()) =1+N(MonteCarlo Average(-( -
(2-13)
However, if N is at all appreciable, 1/S2willnot, differ much from unity so
that this swindle should not significantly improveupon the earlier one unless
the sample size, N, is quite small. Thisagrees with the folklore.-6-
3. Location withGauss/IndependentErrors
Criticalto the swindle in the Gaussian case wastheability to evaluate
(2—li) and (2—12) exactly. These calculations lean heavily on properties of
and S2 in the Gaussian distribution. It is not clear how to successfully
geiera1ize this to arbitrary symmetric unimodal densities, f. However, the
class of distributions given by
y. =uS/v. (3-1)
where u. is unit Gaussian and v. i a positiverandomvariable indeoendent of 1 1
u..,is such that an expression analogous to (2-li) cah be evaluated exactly
and one analogous to (2-12) can be evaluated exactly in some cases and partially
evaluated in all cases. This leads to a swindle that is not as effective as
the one for the Gaussian case, but which is better than simple random sampling.
The family of densities associated with (3-i) is a generalization of the t-
family and contains such meibers as: Cauchy, t, double exponential, logistic
and scale mixtures of Gaussian densities. Conditionally, gven V.y.is
Gaussian with mean zero and variance vT2. We may regard y. as Gaussian with
a randomscale; Andrews and Mallows [1973] give conditions under which a
density has the representation (3.-i).
Letv =(v1,...,vs);then the key idea is thatgiven v,we areback in
muchthe same situation as we were in the pure Gaussian case. The only real
differences are (1) nowthe variances are unequal and (2) wemust eventually
integrate over thedensity of v.Welet
Ev2.y.
.= (v) = (3-2).
Z
1
and
=S2(v)= Zv(y. -(v))2. (33)-7—
Note that (v) and S2(v) can't be computed in real data sincev is not an
observable but in a Monte Carlo study in which v is generated along with
u =(u1,...,uN)Tto produce y, v will be available.
Now instead of knowing the best location-scale invariant estimator of
for the error distribution given by (3-1), we know an even betterestimator,
(v). It is better than the best location—scale invariant estimator because
it uses unobservable information. Thus we will try tocompute the excess
variances ofover that of (v). The formula for this is given in the next
theorem.
Theorem 2: Ifis location and scale invariant estimator, then if the
errors are given y_ (3-1) we have
Var(s) =Var((v))+ E( -(v))2 (34)
=Var((v))+ E((v))2
Proof:
We first show (3-4) and then derive (3-5) from it. Givenv we may compute
=(v)so that we have
((y))2 =((y-))2+ 2 (y -)y+ ()2
and hence taking conditional expectations we get
E[(B(y))2v] =E[((y-))21v]+2 E[ (y -)fv]+ E[()2tvJ
However given v, y and y -yare independent so the middle term vanishes.
Then taking expectations over v we get (3-4). To prove (3—5) we need to show
that
-y(v))2=_______-8--
We have
-y(vfl2v]=E[S2(V)si)Iv] (3-6)
but given v, S(v) is independent of (y —(v))/S(v)so that S2(v) is
independent of Hence (3-6) equals
E[S2(v)Iv] E[ )V))2Iv]
- - (v))2 —
S2(v)
V
Tingexpectations over v proves the result. OED
Tous Theorem 2 to get a swindle we need a little more work, namely we
needto be able to compute Var((vfl. In general this is difficult, but the
following result helps a little.
Theorem 3: Var((v)) E( ). (3-7) (
zv.2
i-I
Proof:Condition on v.
Depending on what the distribution of V. is, the simplification implied
by Theorem 3 may or may not lead to an exact solution. When qv has a chi—
square distributior with q degrees of freedom, then we may obtain an exact
result for Var( (v)). It is:
E(
1
-2 (3-8) q
ii
The swindle flOWcomesin two forms depending on whether or not we have
an exact formula for Var((v)). As before wegiveit for estimating N Var(s).
If Va•r((v)) is known then use:
A A (
Estimatedt (Var(s)) =NVr (v) + N(MonteCarlo Average(. — (v))2)
(3-9)-9-
If Var(v)) cannot be found exactly thenuse
Esttted 1(Var(s)) = N(t1onteCarlo Average(
Zv1
+ tl(oiteCarlo Average( —(v))2) (3-10)
As before, some extra swindlemay be gained from using (3-5) rather than (3-4),
but unless N is small the gains are not likely to beappreciable.
In this case the swindle has two things goingagainst it. 1'ost obviously,
ifE()can't be computed exactly, and must be estimatedby simple random
samplingthen not only are weusing Monte Carlo to estimate the excess variance,
weare alsousing itto estimatea portionofVar((v)). Secondly,even if
Var((v))can becomputedexactly,itmay not beavery largeportionof the
total variance of .Thisis because(v) is a better estimator thanany
locationand scale invariant estimator, i.e.,
Var((v)) <Var(s) (3-11)
forany such estimator . Thusrelativetoany given location—scale invariant
estimatorVarc())may be very small, even relative to the best such estimator.
Because of these problems, the swindle should not beas effective here s it
is for the Gaussian case.-id-
4. The Regression Case
Withthe preparation given inthe previous two sections we may move easily
to the regression case. We first treat the case of Gaussian errors and then
Gauss/Independent errors. The theorems are stated without proof since they
are completely analogous to the corresponding ones for the location se.
Gaussian Errors
If the errors are Gaussian, then we have the following basic result.
Theorai 4: Ifis any rearession and scale invariant estimator and the _____— _______ — ______ ______
'sual least sguares.estimator, then
Cov()Cov(Ls) + E(( -- LS (4-i)
=
Cov(Ls)+ E(()( - - (4-2)
whereS2 is the usual unbiased estimate of2 based on the leastsquares
residual mean-square.
As was true for Theorem 1, (4-1) only requires regression invarianc,
while (4-2) requires both regression and scale invariance.
Since Cov() is given exactly by
Cov(Ls) =(xTxy (4-3)
we are led to the following Monte Carlo s•,indle formula.
Estimated Cov()(XTXY1 +MonteCarlo Average(( - -
(4..4)
As before when Nis small (actually when N-p is small) theremay be sc—e
additionaladvantage tobasingthe swindle on (4-2) rather than (4-1), but (
otherwise•the rnprovement over (4—4) is not likely to be noticeable.—11—
Gauss/Independent Errors
When the errors have the structuregiven by (3-1), wemaydefine, for
each v, these quantities:
(v) (Xv2>x)xT<v2>, (45)
y(v)X(v)
(4-6)
arid
S2(v)(t-pr1 E (y. -y.(v))2 (47)
Where<v2> is the diagonal matrix basedon v2 =(V,...,
Thenwe have the following theorem.
Theorem 5:Ifis y regression and scale invariantestimator andthe
errorsare given (3—1) then
Cov() =Cov((v))+E((-(v))(-(v))T) (4-8)
Cov((v))+E((s)(
-(v))(-(v)) (4-9)
In order to use Theorem 5 to geta swindle we need to be able tocompute
Cov(8(v)). As before, this is generallydifficult, but can be partially
accomplished from the following result.
Theorem 6: Cov((v)) =E(XT<v2>x)_.
(4—10)
Proof:Condition on v.
There do not appear to be toomany cases where E(XT<v2>x)_lcanbe computed
exactly so that either approximations or onte Carloestimates must be used.
Again we get two swindle formulas depending on whatwe use for Cc.vC(v)). If
Cov((v)) is known or can be well approximated, thenuse-12-
Estimated Cov() Cov((v)) +MonteCarlo Averace(( - - (v))T)
(4-11)
If E(X'<V2>X) mustbe estimated byMonte Carlo use
Estimated Cov()= MonteCarloAverage((XT<v2>Xyl)+ —
(4—12)
Monte CarloAverage(( — (v))( — (v)))
Inthe regressioncase,because wewillusually have to use Monte Carlo
is estimate E(XT<v2>X)_l it is likely that the siindlewill not produce much
of an improvement over simple random sampling.
(
(—13-
5. SomeFinal Remarks
Thebasic result that underlies all ofthis is (2-2). This is
aspecial case of a general result that holdsforthe best location invariant
estimatorfor any given distribution (someconditions may be necessary).
This result is given but notproved in the next theorem.
Theorem 7: Ifis any location invariantestimator andis the best
location invariant estimator, then
Var(s)Var(0) +E( - (5-1)
It is evident that if
0.can be ccputed easily, and 'If Var(s) can be
cOmputedexactly then (5—1) providesa basis for a Monte Carloswindle. Un-
fort.uately,neither norits va!'ianc can beeasily Computed for many cases
outside ofthe Gaussian. The use of the3(v) is a compromise for this state
of affairs.
Role of Configurations
In the discussion of this MonteCarlo swindle by Andrews,et al., [1972]
the concept of a Mconfigurationplays a prominent role.In the development
here no such concept arises. Wecomment on this briefly now.
A configuration is a sample (they's) adjusted in a particularway. One
important example of a configuration is
Y-X3 (5-2)
t!e note that aconfiguration is unchancedby the addition of Xy toy for any
choices of y.A scale invariant COnfigUration
is given by
(y- XLJ/S. (53)-14-
The reader is referred to Andre,.s, etal., [1972]for the use of configur-
ations in the derivation of the swindles discussed inthe previous sections.
In our derivation here there is a place forregression or regression and
scale invariant configurations, but theyare not central to the swindle as
such, rather they may be used to make some of thecomDuting more efficient.
In (4—4) we may compute the Monte CarloAverage in one of two ways:
MonteCarlo Average of ((- -
(5...4)
or " ((y-X3)B(y-XEJT) (5-5)
Ifweare going to compare several estimators then froma computational stand
pointitmakes sense to compute
LS'thenform theconfigurations,y -
and then computeon these rather thanon y. This saves a large number of
subtractions. If a regression and scale swindle isgoing to he used then
raLh￿r than
Monte Carlo Average of - -
(5—6)
itis more efficient to form the configurations givenby (5—3) and compute
(5—6) via
-'y- X ,y-
Monte Carlo Average of ((LS S)T) (57)
Similarremarks hold for theconfigurationsthat arise from the swindle in
the non-Gaussian case.

Monte Carlo for Robust Regression: The Swindle Unmasked

Abstract

Similar works

Full text

Available Versions

Research Papers in Economics