Description
For a sample of 601 persons from a magazine survey (Psychology Today)
conducted in July 1969 issue, Fair modeled the 'number of extramarital
affairs' occurence during the last year (zero or some positive value
from 1 to 12) for a proxy of the US population (both males and females).
A tobit model was used at that time to determine the effects of various
demographic variables.
We compute the censored least absolute deviations (Powell, 1984) estimator
using a MIP formulation due to Bilias, Florios and Skouras.
Bilias, Y, Florios, K, and Skouras, S, Exact computation of censored least
absolute deviations estimators. Tech. rep., Athens University of Economics
and Business & National Technical University of Athens, 2013
Fair, R, A theory of extramarital affairs. Journal of Political Economy
86, 45-61, 1978
Powell, JL, Least absolute deviations estimation for the censored
regression model. Journal of Econometrics 25, 303-325, 1984.
Keywords: mixed integer linear programming, statistics, censored regression models,
censored least absolute deviations estimator
Large Model of Type : MIP
Category : GAMS Model library
Main file : clad.gms includes : claddat.gdx
$title Computation of Fair's extramarital Affairs Model Estimates (CLAD,SEQ=397)
$onText
For a sample of 601 persons from a magazine survey (Psychology Today)
conducted in July 1969 issue, Fair modeled the 'number of extramarital
affairs' occurence during the last year (zero or some positive value
from 1 to 12) for a proxy of the US population (both males and females).
A tobit model was used at that time to determine the effects of various
demographic variables.
We compute the censored least absolute deviations (Powell, 1984) estimator
using a MIP formulation due to Bilias, Florios and Skouras.
Bilias, Y, Florios, K, and Skouras, S, Exact computation of censored least
absolute deviations estimators. Tech. rep., Athens University of Economics
and Business & National Technical University of Athens, 2013
Fair, R, A theory of extramarital affairs. Journal of Political Economy
86, 45-61, 1978
Powell, JL, Least absolute deviations estimation for the censored
regression model. Journal of Econometrics 25, 303-325, 1984.
Keywords: mixed integer linear programming, statistics, censored regression models,
censored least absolute deviations estimator
$offText
Set
p 'explanatory variables'
/ RateMar 'How rate marriage, 1-5 scale, 5=very happy ... 1=very unhappy'
Age 'Age, 17.5 to 57'
YearsMa 'No. of Years married, 0.125 to 15'
Intcpt 'intercept' /
T 'sample size (households)' / 1*601 /;
Parameter
X(T,*) 'explanatory and dependent variables'
y(T) 'value of left censored at zero dependent variable';
$gdxIn claddat
$load X
y(T) = X(T,'Affairs');
$if not set normalize_X $set normalize_X 1
$if not set normalize_Y $set normalize_Y 1
* set normalize_X=1/0 and normalize_Y=1/0
* preferable combination is normalize_X=normalize_Y=1
Parameter
delta 'domain for every parameter to be estimated' / 15 /
Xnms(T,p) 'matrix X, normalized all variances equal to 1 if %normalize_X%==1'
mean(p) 'average of X(T.p) over T for mu sigma normalization'
stdev(p) 'stdev of X(T.p) over T etc'
ynms(T) 'vector y, normalized with variance equal to 1 if %normalize_Y%==1'
meanY 'average of y(T) over T for mu sigma nomalization'
stdevY 'stdev of y(T) over T etc'
omega(T) 'tight valid big M coefficient for disjunctive constraints';
mean(p) = sum(T, X(T,p))/card(T);
stdev(p) = sqrt(sum(T, sqr(X(T,p) - mean(p)))/(card(T) - 1));
meanY = sum(T, y(T))/card(T);
stdevY = sqrt(sum(T, sqr(y(T) - meanY))/(card(T) - 1));
Xnms(T,p) = X(T,p);
$if %normalize_X% == 1 Xnms(T,p) = 1;
Xnms(T,p)$stdev(p) = (X(T,p) - mean(p))/stdev(p);
ynms(T) = y(T);
$if %normalize_y% == 1 ynms(T) = 1;
ynms(T)$stdevY = (y(T) - meanY)/stdevY;
omega(T) = delta*sum(p, abs(Xnms(T,p)));
Scalar RHS 'RHS of disjunctive constraints. dependent on %normalise_y%=1 or 0';
RHS = ((0 - meanY)/stdevY)$(%normalize_y% = 1);
Variable
z 'objective function variable'
beta(p) 'beta coefficients'
phi(T) 'substitutes max(X(T,p)*beta(p), 0)'
gamma(T) 'only auxiliary binary variable for phi computation. substitution of max()'
sm(T) 'slack auxiliary variable for obj.fun. computation. substitution of | |'
sp(T) 'surplus auxiliary variable for obj.fun. computation. substitution of | |';
Binary Variable gamma;
Positive Variable sm, sp;
Equation
objfun 'objective function'
con_phi_a(T) 'constraint for phi computation. a'
con_phi_b(T) 'constraint for phi computation. b'
con_phi_c(T) 'constraint for phi computation. c'
con_phi_d(T) 'constraint for phi computation. d'
con_s(T) 'constraint for sm and sp computation';
objfun.. z =e= sum(T,sm(T) + sp(T));
con_phi_a(T).. phi(T) =g= sum(p,beta(p)*Xnms(T,p));
con_phi_b(T).. phi(T) =g= RHS;
con_phi_c(T).. phi(T) =l= sum(p,beta(p)*Xnms(T,p)) + omega(T)*(1 - gamma(T));
con_phi_d(T).. phi(T) =l= RHS + omega(T)*gamma(T);
con_s(T).. ynms(T) - phi(T) + sm(T) - sp(T) =e= 0;
Model CensoredLADPowell84 / all /;
* Depending on the solver box constraints help or do not help
* They are implicitely enforced through omega
* It is preferable to locate the betas in the interior of the box
* Thus just not include the box constraints explicitely
* beta.lo(p) = -delta;
* beta.up(p) = +delta;
option optCr = 0;
solve CensoredLADPowell84 minimizing z using mip;
Parameter
ffbeta(p) 'parameter vector components'
fbeta(p) 'intermediate vector';
ffbeta(p) = beta.l(p);
Alias (p,pp);
if(%normalize_X% = 0,
if(%normalize_Y% = 0,
* 00 case. normalizeX,normalizeY
ffbeta(p) = beta.l(p);
else
* 01 case. normalizeX,normalizeY
fbeta(p) = sum(pp$(stdev(pp) = 0), beta.l(pp)*stdevY + meanY);
fbeta(p)$stdev(p) = (beta.l(p)*stdevY);
ffbeta(p) = fbeta(p);
);
else
if(%normalize_Y% = 0,
* 10 case. normalizeX,normalizeY
fbeta(p) = -sum(pp$stdev(pp), beta.l(pp)*mean(pp)/stdev(pp)) + beta.l(p);
fbeta(p)$stdev(p) = beta.l(p)/stdev(p);
ffbeta(p) = fbeta(p);
else
* 11 case. normalizeX,normalizeY
fbeta(p) = (-sum(pp$stdev(pp), beta.l(pp)*mean(pp)/stdev(pp)) + beta.l(p))*stdevY + meanY;
fbeta(p)$stdev(p) = (beta.l(p)*stdevY)/stdev(p);
ffbeta(p) = fbeta(p);
);
);
display beta.l, ffbeta;
Parameter devRaw(T), zRaw,sRaw;
devRaw(T) = abs(y(T) - max(sum(p, ffbeta(p)*X(T,p)),0));
sRaw = sum(T,devRaw(T));
zRaw = sRaw;
display zRaw;