Code: Select all
'Robust Regression
'Reference: Adapted from the MATLAB code written by James P. LeSage, Dept of Economics, University of Toledo. Original version of the code is available at: http://www.spatial-econometrics.com/regress/robust.m
call robustreg(y,x,4,1)
' ----------------------------------------------------------------------------------------------------
' Program arguments
'-----------------------------------------------------------------------------------------------------
'series y ' dependent variable
'group x ' group of independent variable(s) (including single series)
'scalar model ' 1 = Huber's t function, 2 = Ramsay's E function,
' 3 = Andrew's wave function, 4 = Tukey's biweight
'scalar addc '1 = adds a constant
' ----------------------------------------------------------------------------------------------------
subroutine robustreg(series y,group x,scalar model,scalar addc)
'The following objects will be created through the program.
'If necessary, assign different names to avoid overwriting.
%coef = "beta"
%eq1 = "ols"
%eq2 = "itrw"
%table = "output"
%ser = "weight"
%resid = "resw"
%const = "constant"
'Get the names of dependent and independent variables
%dep = y.@name
%indep = x.@name
'Generate a series for the constant
if addc = 1 then
series {%const} = 1
{%indep}.add {%const}
endif
'Determine the number of independent variables.
!n = {%indep}.@count
'Construct a table to store the results
table(4+!n,5) {%table}
{%table}.setwidth(1:5) 10
{%table}.setlines(a2:e2) +d
{%table}.setlines(a4:e4) +d
{%table}.setformat(D) f.5
{%table}(3,1)= "Variable"
{%table}(3,2)= "Coefficient"
{%table}(3,3)= "t-Statistic"
{%table}(3,4)= "Prob."
{%table}(3,5)= "R-squared"
for !v = 1 to !n
{%table}(4+!v,1) = {%indep}.@seriesname(!v)
next
'Find starting values
equation {%eq1}.ls {%dep} {%indep}
vector {%coef} = @subextract(c,1,1,!n,1)
series {%resid} = resid
!wparm = 2*({%eq1}.@se)^2 'weighting parameter
!scale = @median(@abs({%resid} - @median({%resid})))/.6745 'scale parameter
'Initialize series and control variables
series {%ser} =1
group {%indep}star
vector(!n) {%coef}0
!tol=1
!count=0
'Start the loop
while !tol >.00001
{%coef}0 = {%coef}
{%resid} = {%resid}/!scale
!count=!count+1
if model = 1 then
{%table}(1,1)= "Huber's t function"
{%ser} = !wparm/@abs({%resid})
{%ser} = @recode(@abs({%resid})<=!wparm,1,{%ser})
else if model = 2 then
{%table}(1,1)= "Ramsay's E function"
{%ser} = @exp(-!wparm*@abs({%resid}))
{%ser} = @recode({%resid} =0,1,{%ser})
else if model = 3 then
{%table}(1,1)= "Andrew's wave function"
{%ser} = @sin({%resid}/!wparm)/({%resid}/!wparm)
{%ser} = @recode({%resid}=0,1,@recode({%resid}>@acos(-1)*!wparm,0,{%ser}))
else if model = 4 then
{%table}(1,1)= "Tukey's biweight"
{%ser} = (1-({%resid}/!wparm)^2)^2
{%ser} = @recode({%resid}=0,1,@recode(@abs({%resid})>!wparm,0,{%ser}))
else
{%table}(1,1)= "Tukey's biweight"
{%ser} = (1-({%resid}/!wparm)^2)^2
{%ser} = @recode({%resid}=0,1,@recode(@abs({%resid})>!wparm,0,{%ser}))
endif
endif
endif
endif
'Weighted least squares
series {%dep}star = {%dep}*@sqrt({%ser})
for !i=1 to !n
series {%indep}star!i = {%indep}(!i)*@sqrt({%ser})
{%indep}star.add {%indep}star!i
next
series {%dep}hat=0
equation {%eq2}.ls {%dep}star {%indep}star
{%coef} = @subextract(c,1,1,!n,1)
for !j=1 to !n
{%dep}hat = {%dep}hat+{%coef}(!j)*{%indep}(!j)
next
{%resid} = {%dep}-{%dep}hat
'Recalculate the tolerance level
!tol = @max(@ediv(@abs({%coef}-{%coef}0),@abs({%coef}0)))
wend
'Generate some basic equation output (can be extended further)
stom({%indep}star,{%indep}m)
!evar = @sumsq({%resid})/(@obs({%dep})-!n)
vector {%coef}sig = @sqrt(!evar*@getmaindiagonal(@inverse(@transpose({%indep}m)*{%indep}m)))
vector {%coef}t= @ediv({%coef},{%coef}sig)
'Store the generated output values into the table
for !v=1 to !n
{%table}(4+!v,2) = {%coef}(!v)
{%table}(4+!v,3) = {%coef}t(!v)
{%table}(4+!v,4) = @tdist({%coef}t(!v),@obs({%dep})-!n)
next
{%table}(5,5) = 1-(@sumsq({%resid})/@sumsq({%dep}-@mean({%dep})))
'Final wrap-up
if addc = 1 then
{%indep}.drop {%const}
delete {%const}
endif
delete {%coef}0 {%coef}t {%coef}sig {%eq2} {%indep}m {%indep}star* {%dep}hat {%dep}star
show {%table}
endsub