Predicting Mean and Value-at-risk for a Scale Regression Model
/*--------------------------------------------------------------------------
SAS Sample Library
Name: sevex08.sas
Description: Example Program from SAS/ETS User's Guide,
The SEVERITY Procedure
Title: Predicting Mean and Value-at-risk for a Scale Regression Model
Product: SAS/ETS Software
Keys: Severity Distribution Modeling
PROC: SEVERITY
Notes:
----------------------------------------------------------------------------*/
/*--------- Define distribution functions that compute the mean ----------*/
proc fcmp library=sashelp.svrtdist outlib=work.means.scalemod;
function BURR_MEAN(x, Theta, Alpha, Gamma);
if not(Alpha * Gamma > 1) then
return (.); /* first moment does not exist */
return (Theta*gamma(1 + 1/Gamma)*gamma(Alpha - 1/Gamma)/gamma(Alpha));
endsub;
function EXP_MEAN(x, Theta);
return (Theta);
endsub;
function GAMMA_MEAN(x, Theta, Alpha);
return (Theta*Alpha);
endsub;
function GPD_MEAN(x, Theta, Xi);
if not(Xi < 1) then
return (.); /* first moment does not exist */
return (Theta/(1 - Xi));
endsub;
function IGAUSS_MEAN(x, Theta, Alpha);
return (Theta);
endsub;
function LOGN_MEAN(x, Mu, Sigma);
return (exp(Mu + Sigma*Sigma/2.0));
endsub;
function PARETO_MEAN(x, Theta, Alpha);
if not(Alpha > 1) then
return (.); /* first moment does not exist */
return (Theta/(Alpha - 1));
endsub;
function STWEEDIE_MEAN(x, Theta, Lambda, P);
return (Theta* Lambda * (2 - P) / (P - 1));
endsub;
function TWEEDIE_MEAN(x, P, Mu, Phi);
return (Mu);
endsub;
function WEIBULL_MEAN(x, Theta, Tau);
return (Theta*gamma(1 + 1/Tau));
endsub;
quit;
/*-------------- Simulate a lognormal sample -------------*/
data test_sev8(keep=y x1-x5
label='A Lognormal Sample Affected by Regressors');
array x{*} x1-x5;
array b{6} _TEMPORARY_ (1 0.75 3 -1 0.25 5);
call streaminit(45678);
label y='Response Influenced by Regressors';
Sigma = 0.25;
do n = 1 to 500;
Mu = b(1); /* log of base value of scale */
do i = 1 to dim(x);
x(i) = rand('NORMAL');
Mu = Mu + b(i+1) * x(i);
end;
y = exp(Mu) * rand('LOGNORMAL')**Sigma;
output;
end;
run;
/*----- Fit all distributions and generate scoring functions ------*/
options cmplib=work.means;
proc severity data=test_sev8 outest=est print=all plots=none;
loss y;
scalemodel x1-x5;
dist _predefined_ stweedie;
outscorelib outlib=scorefuncs commonpackage;
run;
/*--- Simulate scoring data ---*/
data reginput(keep=x1-x5);
array x{*} x1-x5;
call streaminit(9041);
do n=1 to 10;
do i = 1 to dim(x);
x(i) = rand('NORMAL');
end;
output;
end;
run;
/*--- Set VaR level ---*/
%let varLevel=0.975;
/*--- Compute scores (mean and var) for the ---
--- scoring data by using the scoring functions ---*/
data scores;
array x{*} x1-x5;
set reginput;
igauss_mean = sev_mean_igauss(., x);
igauss_var = sev_quantile_igauss(&varLevel, x);
logn_mean = sev_mean_logn(., x);
logn_var = sev_quantile_logn(&varLevel, x);
run;
/*--- Compute scores (mean and var) for the ---
--- scoring data by using the OUTEST= data set ---*/
data scoresWithOutest(keep=x1-x5 igauss_mean igauss_var logn_mean logn_var);
array _x_{*} x1-x5;
array _xparmIgauss_{5} _temporary_;
array _xparmLogn_{5} _temporary_;
retain _Theta0_ Alpha0;
retain _Mu0_ Sigma0;
*--- read parameter estimates for igauss and logn models ---*;
if (_n_ = 1) then do;
set est(where=(upcase(_MODEL_)='IGAUSS' and _TYPE_='EST'));
_Theta0_ = Theta; Alpha0 = Alpha;
do _i_=1 to dim(_x_);
if (_x_(_i_) = .R) then _xparmIgauss_(_i_) = 0;
else _xparmIgauss_(_i_) = _x_(_i_);
end;
set est(where=(upcase(_MODEL_)='LOGN' and _TYPE_='EST'));
_Mu0_ = Mu; Sigma0 = Sigma;
do _i_=1 to dim(_x_);
if (_x_(_i_) = .R) then _xparmLogn_(_i_) = 0;
else _xparmLogn_(_i_) = _x_(_i_);
end;
end;
set reginput;
*--- predict mean and VaR for inverse Gaussian ---*;
* first compute X'*beta for inverse Gaussian *;
_xbeta_ = 0.0;
do _i_ = 1 to dim(_x_);
_xbeta_ = _xbeta_ + _xparmIgauss_(_i_) * _x_(_i_);
end;
* now compute scale for inverse Gaussian *;
_SCALE_ = _Theta0_ * exp(_xbeta_);
igauss_mean = igauss_mean(., _SCALE_, Alpha0);
igauss_var = igauss_quantile(&varLevel, _SCALE_, Alpha0);
*--- predict mean and VaR for lognormal ---*;
* first compute X'*beta for lognormal*;
_xbeta_ = 0.0;
do _i_ = 1 to dim(_x_);
_xbeta_ = _xbeta_ + _xparmLogn_(_i_) * _x_(_i_);
end;
* now compute Mu=log(scale) for lognormal *;
_MU_ = _Mu0_ + _xbeta_;
logn_mean = logn_mean(., _MU_, Sigma0);
logn_var = logn_quantile(&varLevel, _MU_, Sigma0);
run;
proc compare data=scoresWithOutest compare=scores
crit=1.0e-12;
run;