Proc IML and Call NLPNRR/NLPNRA for Maximum Likelihod

thaisjensen — Wed, 30 May 2012 16:44:15 GMT

Hey SAS-community

I have a problem with maximizing a likelihood function within proc iml using call nlpnra and nlpnrr. The data I am looking at is event/survival data where I am modeling defaults and at each event time will have to take into account observations in the risk set. This is what is specified in the function statement below.

The program outputs a value for the objective function(l), for the initial values of the beta-vector(betaT) but does however not try to change it for a better fit. I have tried running nlpnra and nlpnrr calls, but get two different warning/error messages.

Below a sketch of the code I am using is included with a simple data example. The problem is however the same when I run it on my larger data set. Is impossible to get at least some rough estimates without specifying both the gradient and the hessian?

Thank you very much for your help.

-- ERRORS --

NLPNRR
ERROR: ITEM \_pages_\_page_#5 does not exist.
NOTE: XCONV convergence criterion satisfied.
NOTE: At least one element of the (projected) gradient is greater than 1e-3.

NLPNRA

ERROR: NEWRAP Optimization cannot be completed.

WARNING: Optimization routine cannot improve the function value.

-- CODE --

dm "cle out";

dm "cle log";

proc iml;

/* ---------------------------------------------------------------------------------------------------------------------*/

/* Simple example: 9 observations, 2 covaraites, 6 of the 9 firms defauklts */

x={1 1, 1 2, 1 3, 1 4, 1 5, 1 6, 1 7, 1 8, 1 9}; /* 9x2 vector - design matrix 2 covaraites where first i constant*/

d={0, 0.4, 3, 4, 6, 7}; /* 6x1 vector - default times of defaulting firms */

x_d={1 1, 1 3, 1 5, 1 6, 1 7, 1 8}; /* 6*2 vector - covariates of defaulting firms*/

time={0, 0.2, 0.4, 0.7, 3, 4, 6, 7, 8}; /* 9x1 vector - observation tim*/

default={1, 0, 1, 0, 1, 1, 1, 1, 0}; /* 9x1 vector - default at observation time*/

/* ---------------------------------------------------------------------------------------------------------------------*/

p = ncol(x); /* Count the number of parameters to be included */

n = nrow(x); /* count the number of observations */

dn = nrow(d); /* counts the number of distinct defaults */

l=0; /* initializes the loglikelihood OUTSIDE ALL LOOPS */

S_0=j(1,1,0); /* initiates S_0 - aalen p 149 */

/* -- can manually be changed ----- */

betaT={0 0}; /* initializes betaT row vectors */

/* ---------------------------------------------------------------------------------------------------------------------*/

start f_loglike(betaT) global(n,dn,x,time,d,x_d,l);

beta=t(betaT); /* transposes the initalized vector - becomes a COLUMN vector */

/* enters a loop from first default time to last default time - 6 loops in total */

do j=1 to dn;

S_0=0; /* resets S back 0 for next default runthrough: j=1 to dj loop*/

dj = d; /* stores jth default time into dj */

dj_round = ceil(dj+0.000001); /* rounds up to nearest interger, where quarter ends */

/* IDEA: create a vector that is equal 1 if observation belongs to risk set for the particular default */

/* if it belongs to risk set then its end_global_frac time is between di and di_round */

dj_vector=j(n,1,dj); /*WITHIN Default-LOOP. Reinitilalizes each time loops runs through*/

dj_round_vector=j(n,1,dj_round); /*WITHIN Default-LOOP. Reinitilalizes each time loops runs through*/

t1= dj_vector <= time; /* Vector equalling one when observation when end_global_frac os greater then or equal to jth default time DIMENSIONS: nx1/*/

t2= time <= dj_round_vector; /* Vector equalling one when observation when end_global_frac os less than or equal to jth default time DIMENSIONS: nx1/*/

y= t1#t2; /* Vector of 1s for firms at risk at the particular jth default time and otherwize zero DIMENSIONS: nx1 */

do i=1 to n; /* summing over all n - only positive values for firms at risk*/

x_i=x[i,]; /* taking out the covariates for the ith firms if at risk 0 */

si=exp(x_i*beta); /* multiplies with beta vector and set in exponents */

si2=y#si; /* muliplies with ith element of t3 vector which is 1 if beloning to the risk set and zero otherwise */

S_0=S_0+si; /* sums over the risk set*/

end;

x_dj=x_d[j,]; /* stores j th defaulting firms covariates at default time into row vector */

xjbeta=x_dj*beta; /* multiplies row vector with column vector to form jth firms contrubution */

iloglike =xjbeta-log(S_0); /* iloglikelihood for jth default substracted S at time t_j default time */

l = l + iloglike; /* adds each seperate log likelihood elements to the agg. likelihood to be maximized*/

end;

return(l);

finish f_loglike;

/* ---------------------------------------------------------------------------------------------------------------------*/

/* http://support.sas.com/documentation/cdl/en/imlug/59656/HTML/default/viewer.htm#nonlinearoptexpls_sect10.htm */

optn = {1 4};

call nlpnra(xr,rc,"f_loglike",betaT,optn);

quit;

Re: Proc IML and Call NLPNRR/NLPNRA for Maximum Likelihod

Rick_SAS — Wed, 30 May 2012 18:28:36 GMT

To optimize a function, it must have a local maximum. Your function is unbounded.

Re: Proc IML and Call NLPNRR/NLPNRA for Maximum Likelihod

thaisjensen — Thu, 31 May 2012 10:43:11 GMT

Thanks Rick for the swift reply. I see why that will cause the optimizer to not work. I'll work at getting the right likelihood specified.

topic Proc IML and Call NLPNRR/NLPNRA for Maximum Likelihod in SAS/IML Software and Matrix Computations

Proc IML and Call NLPNRR/NLPNRA for Maximum Likelihod

Re: Proc IML and Call NLPNRR/NLPNRA for Maximum Likelihod

Re: Proc IML and Call NLPNRR/NLPNRA for Maximum Likelihod