@ Space allocation for program. Do not change. @ new , 30000; @ Output to screen? (in addition to file) @ screen on; @ Name of the output file, reset=overwrite existing file, on=add to it @ output file=c:\temp2\output.txt reset; @ Put a title for your output in the quotes below. @ print "Mix logit with 1 fixed, 1 normal, 2 uniform, and 2 triangular."; print ; @ 1 to check inputs for conformability and consistency, else 0. @ VERBOSE = 1; @ Number of people. (integer) @ NP = 84; @ Number of observations. (integer) @ NOBS = 1004; @ Maximium number of alternatives. (integer) @ NALT = 4; @ Load or create XMAT (contain all of the explanatory variables and the@ @ censoring variable (if needed), YVEC (the dependent variable), @ @ and TIMES (giving the number of choice situations for each person.)@ @----------------------------------------------------------------@ load XMAT[NOBS,28] = c:\temp2\xmat.asc; load YVEC[NOBS,1] = c:\temp2\yvec.asc; load TIMES[NP,1] = c:\temp2\times.asc; @ Variables are now: price, contract length, LL, KK+KL, UL, TOD, seasonal.@ @----------------------------------------------------------------@ @ Number of variables in XMAT. (integer) @ NVAR = 7; @ Give the number of explanatory variables that have fixed coefficients. @ @ (integer) @ NFC = 1; @ NFCx1 vector to identify the variables in XMAT that have fixed @ @ coefficients. (column vector) @ IDFC = { 1 }; @ Give the number of explanatory variables that have normally @ @ distributed coefficients. (integer) @ NNC = 1; @ NNCx1 vector to identify the variables in XMAT that have normally @ @ distributed coefficients. (column vector) @ IDNC = { 2 }; @ Give the umber of explanatory variables that have uniformly @ @ distributed coefficients @ NUC = 2; @ NUCx1 vector to identify the variables in XMAT that have uniformly @ @ distributed coefficients. (column vector) @ IDUC = { 3, 4 }; @ Give the number of explanatory variables that have triangularly @ @ distributed coefficients @ NTC = 2; @ NTCx1 vector to identify the variables in XMAT that have triangularly @ @ distributed coefficients. (column vector) @ IDTC = { 6, 7 }; @ Give the number of explanatory variables that have log-normally @ @ distributed coefficients. (integer) @ NLC = 0; @ NLCx1 vector to identify the variables in XMAT that have @ @ log-normally distributed coefficients. (column vector) @ IDLC = { 0 }; @ 1 if all people do not face all NALT alternatives, else 0. @ CENSOR = 0; @ Identify the variable in XMAT that identifies the alternatives @ @ each person faces. (integer) @ IDCENSOR = 0; @ If you want to weight the observation, set DOWT to 1, else 0 @ DOWT = 0; @ If DOWT=1, load or create an NPx1 vector of weights. @ WEIGHT=0; @ 1 to print out the diagonal elements of the Hessian, else 0. @ PRTHESS = 0; @ 1 to rescale any variable in XMAT, else 0. @ RESCALE = 0; @ If RESCALE = 1 create a q x 2 matrix to id which variables to @ @ rescale and by what multiplicative factor, else make RESCLMAT @ @ a dummy value. @ RSCLMAT = { 0 }; @ Starting values for the parameters. @ @ (NFC+(2*NNC)+(2*NUC)+(2*NTC)+(2*NLC)) x 1 vector of starting values. @ @ Order is NFC fixed coefficients, followed by mean and standard deviation of each @ @ of NNC normal coefficients, followed by mean and standard deviation of each of NUC@ @ uniformly distributed coefficients, etc for the triangularly and log-normally @ @ distributed coefficients.@ B = { -.90, -.21, .39, 2.40, 2.94, 1.55, 1.99, -8.68, 4.59, -8.73, 4.92 }; @ 1 to constrain any parameter to its starting value, else 0. @ CZVAR = 0; @ If CZVAR = 1, create a vector of zeros and ones of the same @ @ dimension as B identifying which parameters are to vary, else make @ @ BACTIVE = { 0 }. (column vector) @ BACTIVE = { 0 }; @ 1 to constrain any of the error components to equal each other, @ @ else 0. (integer) @ CEVAR = 0; @ If CEVAR=1, create a matrix of ones and zeros to constrain the @ @ necessary random parameters to be equal, else make EQMAT=0. (matrix)@ EQMAT = { 0 }; @ Number of repetitions. NREP must be >= 1. (integer) @ NREP = 125; @ Choose random or halton draws. DRAWS=1 for random, 2 for Halton.@ DRAWS=2; @ If DRAWS=1, set seed to use in random number generator. SEED1 must be >= 0. (integer) @ SEED1 = 46; @ If DRAWS=2, specify the following.@ @ HALT=1 to create Halton draws, HALT=2 to use previously created Halton draws. @ HALT=1; @ If HALT=1, set SAVH=1 if you want to save the Halton draws that are created; 0 otherwise. @ SAVH=0; @ HMNAME = path and name for file of halton sequences. Be sure to put double \\ where@ @ single \ appears in the path. @ @ If HALT=1 and SAVH=1, the Halton draws that are created are saved to this file.@ @ If HALT=2, the Haton draws that were previously saved to this file are loaded.@ @ Note that, if HALT=2, the sequences must meet the specs of the current model,@ @ that is, the same NREP, NOBS, number of random coefficients, and distribution for each coefficient.@ HMNAME="c:\\temp2\\hm125.fmt"; @ Maximum number of iterations in the maximization. (integer) @ NITER = 100; @ Tolerance for convergence. (small decimal) @ EPS = 1.e-4; @ Identify the optimization routine: 1 Paul Ruud's routine @ @ 2 for Gauss's maxlik @ OPTIM = 1; @ Specify the optimization algorithm/hessian calculation.(integer) @ @ Options with OPTIM=1: 1 for bhhh, 2 for numerical (newton-raphson) @ @ Options with OPTIM=2: 1 for steepest descent, 2 for bfgs, 3 for dfp @ @ 4 for numerical (newton-raphson), 5 for bhhh, 6 for prcg. @ METHOD = 1; @ 1 if OPTIM=1 and want to use modified iteration procedure; else 0. @ MNR = 1; @ If OPTIM=1, set STEP to the step length the maximization routine @ @ should initially try. @ STEP = 1; @ If OPTIM=1, then set ROBUST=1 if you want robust standard errors, or@ @ ROBUST=0 if you want regular standard errors. With OPTIM=2, this option doesn't exist.@ ROBUST=1; @======================================================================@ @@@ You should not need to change anything below this line @@@ @======================================================================@ @ Create global for the number of estimated variables @ NEVAR = NFC+(NNC+NUC+NTC+NLC)*2; @ Check inputs if VERBOSE=1 @ if ((VERBOSE /= 1) and (VERBOSE /= 0)); print "VERBOSE must 0 or 1."; print "Program terminated."; stop; endif; if VERBOSE; pcheck; else; print "The inputs have not been checked since VERBOSE=0."; print; endif; @ Rescale the variables. @ if RESCALE; j = rows(RSCLMAT); i = 1; if VERBOSE; if (RSCLMAT[.,1] > NVAR); print "RSCLMAT identifies a variable that is not in the data set."; print "Program terminated."; stop; endif; print "Rescaling Data:"; print " Variable Mult. Factor"; do while i <= j; RSCLMAT[i,1] RSCLMAT[i,2]; XMAT[.,(RSCLMAT[i,1]-1)*NALT+1:RSCLMAT[i,1]*NALT] = XMAT[.,(RSCLMAT[i,1]-1)*NALT+1:RSCLMAT[i,1]*NALT] * RSCLMAT[i,2]; i = i + 1; endo; print " "; else; do while i <= j; XMAT[.,(RSCLMAT[i,1]-1)*NALT+1:RSCLMAT[i,1]*NALT] = XMAT[.,(RSCLMAT[i,1]-1)*NALT+1:RSCLMAT[i,1]*NALT] * RSCLMAT[i,2]; i = i + 1; endo; endif; endif; @ Print out starting values if VERBOSE = 1. @ if VERBOSE; print "There are " NP " respondents and " NOBS " observations."; print "The model has" NFC " fixed coefficients, for variables" IDFC; print " and" NNC " normally distributed coefficients, for variables" IDNC; print " and" NUC " uniformly distributed coefficients, for variables" IDUC; print " and" NTC " triangularly distributed coefficients, for variables" IDTC; print " and" NLC " log-normally distributed coefficients, for variables" IDLC; print; if CZVAR; print "Starting values:" ; print B; print ; print "Parameters that are estimated (1) or held at starting value (0):"; print BACTIVE; print ; endif; if (CZVAR == 0); print "Starting values:"; print B; print; print "All parameters are estimated; none is held at its starting value."; print; endif; endif; @ Create new B and BACTIVE with reduced lengths if @ @ user specifies equality constraints. @ if CEVAR; if CZVAR; BACTIVE = EQMAT * BACTIVE; BACTIVE = BACTIVE .> 0; endif; B = (EQMAT * B) ./ (EQMAT * ones((NFC+2*NNC+2*NUC+2*NTC+2*NLC),1)); if VERBOSE and CZVAR; print "Some parameters are constrained to be the same."; print "Starting values of constrained parameters:"; print B; print ; print "Constrained parameters that are estimated (1) or held at starting value (0):"; print BACTIVE; print; endif; if VERBOSE and (CZVAR == 0); print "Some parameters are constrained to be the same."; print "Starting values of constrained parameters:"; print B; print; print "All constrained parameters are estimated; none is held at its starting value."; print; endif; endif; @ Describing random terms. @ if VERBOSE; if DRAWS == 1; print "Random draws are used."; print "Random error terms are based on:"; print "Seed: " SEED1; print; elseif DRAWS == 2; print "Halton draws are used."; print; else; print "DRAWS must be 1 or 2. Program terminated."; stop; endif; print "Repetitions: " NREP; endif; /* Number of random coefficients or 1, whichever is higher. */ NECOL = maxc( ones(1,1) | (NNC+NUC+NTC+NLC) ); if DRAWS == 2; /* CREATE OR LOAD HALTON SEQUENCES. */ if HALT == 2; loadm hm = ^HMNAME; print "Loaded Halton sequences."; elseif HALT == 1; @ Create the Halton sequence @ print "Creating Halton sequences ...."; /* Provide prime number vector */ prim = { 2 3 5 7 11 13 17 19 23 29 31 37 41 43 47 53 59 61 67 71 73 79 83 89 97 101 103 107 109 113 }; print "Halton sequences are based in primes: " prim[1,1:NECOL]; print; /* Develop a halton matrix (hm) of size nrep x (nobs*necol); The first 'nobs' columns provide halton numbers for first dimension of integration, the second 'nobs' columns provide the halton numbers for the second dimension, and so on; cdfinvn is the inverse cumulative standard normal distribution function; PROC halton creates sequence; PROC cdfinvn takes inverse normal: windows version of gauss has this as a command. For non-normal distributions, use inverse of other distribution. */ h = 1; hm = {}; do while h <= (NECOL); hm1 = halton(10+nrep*np,prim[h]); if h <= NNC or h > (NNC+NUC+NTC); @ Normal and lognormal distributions. @ @ If there are no random coefs, st NECOL=1 while NNC=NUC=NTC=NLC=0, @ @ then, hm1 will be treated as normally distributed.@ hm1 = cdfinvn(hm1); @ The inverse-normal proc produces very extreme values sometimes. This truncates.@ hm1=hm1.*(hm1 .le 10) + 10 .* (hm1 .gt 10); hm1=hm1.*(hm1 .ge -10) -10 .* (hm1 .lt -10); endif; if h > NNC and h <= (NNC+NUC); @Uniform distribution.@ hm1=hm1.*2 - 1; endif; if h > (NNC+NUC) and h <= (NNC+NUC+NTC); @Triangular distribution.@ hm1= (sqrt(2.*hm1)-1).*(hm1 .<= 0.5) + (1-sqrt(2.*(1-hm1))) .* (hm1 .> 0.5); endif; hm = hm~hm1[11:rows(hm1),1]; h=h+1; endo; print "Finished Halton sequences."; if SAVH; save ^HMNAME = hm; print "Saved Halton sequences."; endif; else; print "HALT must equal 1 or 2. Job terminated."; stop; endif; print "Number of rows in Halton matrix: " rows(hm); print "Number of columns in Halton matrix: " cols(hm); endif; @ END OF HALTON DRAWS @ @ IDA is a vector that identifies the variables with normal, uniform or triangular coefficients.@ @ It is used in LL, gradient and forecst procs. @ IDA={0}; if NNC>0; IDA=IDNC; endif; if NNC==0 and NUC>0; IDA=IDUC; endif; if NNC>0 and NUC>0; IDA=IDA|IDUC; endif; if (NNC+NUC)==0 and NTC>0; IDA=IDTC; endif; if (NNC+NUC)>0 and NTC>0; IDA=IDA|IDTC; endif; @ Create dependent variable permutation matrix for the gradient @ YPERM=zeros(NOBS,NALT); i = 1; do while i<=NOBS; YPERM[i,YVEC[i,1]] = 1; i = i + 1; endo; @ Initialize the STEP to twice its value for Paul Ruud's routine. @ if (OPTIM==1); STEP = STEP*2; endif; @ Set up and do maximization @ if OPTIM == 1; beta = domax(&ll,&gr,B,BACTIVE); print; print "Remember: Signs of standard deviations are irrelevant."; print "Interpret them as being positive."; print; endif; if OPTIM == 2; library maxlik,pgraph; #include maxlik.ext; maxset; _max_GradTol = EPS; _max_GradProc = &gr; _max_MaxIters = NITER; _max_Algorithm = METHOD; if DOWT; __weight = WEIGHT; endif; if CZVAR; _max_Active = BACTIVE; endif; {beta,f,g,cov,ret} = maxlik(XMAT,0,&ll,B); call maxprt(beta,f,g,cov,ret); print; print "Remember: Signs of standard deviations are irrelevant."; print "Interpret them as being positive."; print; if (CZVAR == 0); print "gradient(hessian-inverse)gradient is:" ((g/_max_FinalHess)'g); else; g = selif(g,BACTIVE); print "gradient(hessian-inverse)gradient is:" ((g/_max_FinalHess)'g); endif; if PRTHESS; print "diagonal of hessian:" ((diag(_max_FinalHess))'); print; endif; endif; @ THIS IS THE END OF THE MAIN PROGRAM. @ @ PROCS ll, gr, gpnr, expand, domax, pcheck follow.@ /* LOG-LIKELIHOOD FUNCTION */ proc ll(b,x); @ Relies on the globals: NALT, NOBS, NP, NREP, NECOL, SEED1, CENSOR @ @ NFC, IDFC, NNC, NUC,NTC, IDA, NLC, IDLC @ @ CZVAR, BACTIVE, CEVAR, EQMAT, XMAT, DRAWS @ local c, k, n, t, km, kmm, rd, seed2; local v, ev, p0, p00, err, mm; if CEVAR; b = EQMAT' * b; @ Expand b to its original size. @ endif; v = zeros(NOBS,NALT); @ Argument to logit formula @ p0 = zeros(NP,1); @ Simulated probability @ if CENSOR; c = (IDCENSOR-1)*NALT; @ Location of censor variable @ endif; k = 1; do while k <= NFC; @ Adds variables with fixed coefficients @ km = (IDFC[k,1]-1)*NALT; v = v + b[k] .* XMAT[.,(km+1):(km+NALT)]; k = k+1; endo; seed2 = SEED1; rndseed seed2; rd = 0; n = 1; do while n <= NP; if DRAWS == 1; err = rndn(NREP,NECOL); if NUC > 0; err[.,(NNC+1):(NNC+NUC)] = (rndu(NREP,NUC) .* 2) -1; endif; if NTC > 0; mm = rndu(NREP,NTC) ; err[.,(NNC+NUC+1):(NNC+NUC+NTC)] = (sqrt(2.*mm)-1).*(mm .<= 0.5) + (1-sqrt(2.*(1-mm))) .* (mm .> 0.5); endif; endif; if DRAWS == 2; err=hm[(NREP*(n-1)+1):(NREP*(n-1)+NREP),.] ; endif; @ err has NREP rows and NECOL columns for each person. @ p00 = ones(NREP,1); t=1; do while (t<=TIMES[n]); kmm = rd + t; ev = v[kmm,.]; k = 1; do while k <= NNC+NUC+NTC; @ Adds variables with normal, uniform and tringular coefficients @ km = (IDA[k,1]-1)*NALT; ev = ev + (b[NFC+(2*k)-1] + (b[NFC+(2*k)] .* err[.,k])) .* XMAT[kmm,(km+1):(km+NALT)]; k = k+1; endo; k = 1; do while k <= NLC; @ Adds variables with log-normal coefficients @ km = (IDLC[k,1]-1)*NALT; ev = ev + exp(b[NFC+(2*(NNC+NUC+NTC))+(2*k)-1] + (b[NFC+(2*(NNC+NUC+NTC))+(2*k)] .* err[.,(NNC+NUC+NTC+k)])) .* XMAT[kmm,(km+1):(km+NALT)]; k = k+1; endo; ev = exp(ev); if CENSOR; p00 = p00.*ev[.,YVEC[kmm,1]]./sumc((ev .* XMAT[kmm,(c+1):(c+NALT)])'); else; p00 = p00.*ev[.,YVEC[kmm,1]]./(sumc(ev')); endif; t = t+1; endo; rd = rd + TIMES[n]; p0[n,1] = meanc(p00); n = n + 1; endo; retp(ln(p0)); endp; /* GRADIENT PROCEDURE */ proc gr(b,x); @ Relies on the globals: NALT, NOBS, NP, NREP, NECOL, SEED1, CENSOR @ @ NFC, IDFC, NNC, NUC, NTC, IDA, NLC, IDLC @ @ CZVAR, BACTIVE, CEVAR, EQMAT, XMAT, YPERM, DRAWS @ local c, g, k, n, t, km, kmm, rd, seed2, y; local denom, der, v, ev, p0, p00, p1, err, mm; if CEVAR; b = EQMAT' * b; @ Expand b to its original size. @ endif; if CENSOR; c = (IDCENSOR-1)*NALT; @ Location of censor variable @ endif; v = zeros(NOBS,NALT); @ Argument to logit formula @ p0 = zeros(NP,1); @ Simulated probability @ der = zeros(NP,NEVAR); @ Jacobian matrix @ k = 1; do while k <= NFC; @ Adds variables with fixed coefficients @ km = (IDFC[k,1]-1)*NALT; v = v + b[k] .* XMAT[.,(km+1):(km+NALT)]; k = k+1; endo; seed2 = SEED1; rndseed seed2; rd = 0; n = 1; do while n <= NP; if DRAWS == 1; err = rndn(NREP,NECOL); if NUC > 0; err[.,(NNC+1):(NNC+NUC)] = (rndu(NREP,NUC) .* 2) -1; endif; if NTC > 0; mm = rndu(NREP,NTC) ; err[.,(NNC+NUC+1):(NNC+NUC+NTC)] = (sqrt(2.*mm)-1).*(mm .<= 0.5) + (1-sqrt(2.*(1-mm))) .* (mm .> 0.5); endif; endif; if DRAWS == 2; err=hm[(NREP*(n-1)+1):(NREP*(n-1)+NREP),.] ; endif; @ err has NREP rows and NECOL columns for each person. @ p00 = ones(NREP,1); g = zeros(NREP,NEVAR); t=1; do while (t<=TIMES[n]); kmm = rd + t; ev = v[kmm,.]; k = 1; do while k <= NNC+NUC+NTC; @ Adds variables with normal,uniform, and triangular coefficients @ km = (IDA[k,1]-1)*NALT; ev = ev + (b[NFC+(2*k)-1] + (b[NFC+(2*k)] .* err[.,k])) .* XMAT[kmm,(km+1):(km+NALT)]; k = k+1; endo; k = 1; do while k <= NLC; @ Adds variables with log-normal coefficients @ km = (IDLC[k,1]-1)*NALT; ev = ev + exp(b[NFC+(2*(NNC+NUC+NTC))+(2*k)-1] + (b[NFC+(2*(NNC+NUC+NTC))+(2*k)] .* err[.,(NNC+NUC+NTC+k)])) .* XMAT[kmm,(km+1):(km+NALT)]; k = k+1; endo; ev = exp(ev); if CENSOR; denom = sumc((ev .* XMAT[kmm,(c+1):(c+NALT)])'); p00 = p00 .* ev[.,YVEC[kmm,1]]./denom; p1 = (ev.*XMAT[kmm,(c+1):(c+NALT)])./denom; else; denom = sumc(ev'); p00 = p00 .* ev[.,YVEC[kmm,1]]./denom; p1 = ev./denom; endif; y = YPERM[kmm,.]-p1; k = 1; do while k<=NFC; km = (IDFC[k,1]-1)*NALT; g[.,k] = g[.,k] + sumc((XMAT[kmm,km+1:km+NALT].*y)'); k = k + 1; endo; k = 1; do while k<=NNC+NUC+NTC; km = (IDA[k,1]-1)*NALT; g[.,NFC+(2*k)-1] = g[.,NFC+(2*k)-1] + sumc((XMAT[kmm,km+1:km+NALT].*y)'); g[.,NFC+(2*k)] = g[.,NFC+(2*k)] + sumc((err[.,k].*XMAT[kmm,km+1:km+NALT].*y)'); k = k + 1; endo; k = 1; do while k<=NLC; km = (IDLC[k,1]-1)*NALT; g[.,NFC+2*(NNC+NUC+NTC)+(2*k)-1] = g[.,NFC+2*(NNC+NUC+NTC)+(2*k)-1] + sumc((XMAT[kmm,km+1:km+NALT].*y)'); g[.,NFC+2*(NNC+NUC+NTC)+(2*k)] = g[.,NFC+2*(NNC+NUC+NTC)+(2*k)] + sumc((err[.,NNC+NUC+NTC+k].*XMAT[kmm,km+1:km+NALT].*y)'); k = k + 1; endo; t = t+1; endo; km = NFC+2*(NNC+NUC+NTC); @ multiply correction term for log-normals @ k = 1; do while k<=NLC; g[.,km+(2*k)-1:km+(2*k)] = g[.,km+(2*k)-1:km+(2*k)] .* exp(b[km+(2*k)-1,1]+b[km+(2*k),1].*err[.,NNC+NUC+NTC+k]); k = k + 1; endo; p0[n,1] = meanc(p00); der[n,.] = (meanc(p00.*g))'; rd = rd + TIMES[n]; n = n + 1; endo; if CEVAR; der = (der./p0) * EQMAT'; else; der = der./p0; endif; retp(der); endp; /* GRADIENT FOR PAUL RUUD'S ROUTINE WHEN USING NEWTON-RAPHSON*/ /* USE WHEN OPTIM == 1 AND METHOD == 2 */ proc gpnr(b); @ Relies on global: XMAT @ local grad; if DOWT; grad = WEIGHT .* gr(b,XMAT); else; grad = gr(b,XMAT); endif; retp(sumc(grad)); endp; /* EXPANDS THE DIRECTION VECTOR; ALLOWS PARAMETERS TO STAY AT STARTING */ /* VALUES; HELPER PROCEDURE FOR &domax */ proc expand( x, e ); local i,j; i = 0; j = 0; do while i < rows(e); i = i + 1; if e[i]; j = j + 1; e[i] = x[j]; endif; endo; if j/=rows(x); "Error in expand."; stop; endif; retp( e ); endp; /* MAXIMIZATION ROUTINE COURTESY OF PAUL RUUD */ proc domax( &f, &g, b, bactive ); @ Relies on the globals: CZVAR, EPS, METHOD, NITER, NOBS, @ @ PRTHESS, XMAT, NP, NEVAR @ local f:proc, g:proc; local direc, first, grad, sgrad, hesh, ihesh, lambda; local nb, printer, repeat, step1, wtsq; local _biter, _f0, _f1, _fold, _tol; _tol = 1; _biter = 0; nb = seqa(1,1,NEVAR); _f0 = f( b, XMAT ); if DOWT; _f0 = sumc(WEIGHT.*_f0); wtsq = WEIGHT .^ (.5); else; _f0 = sumc(_f0); endif; format /m1 /rdt 16,8; print; print; print; do while (_tol > EPS or _tol < 0) and (_biter < NITER); _biter = _biter + 1; print "=========================================================================="; print " Iteration: " _biter; print " Function: " _f0; grad = g( b, XMAT ); if (METHOD == 1); if DOWT; grad = wtsq.*grad; hesh = grad'grad; grad = wtsq.*grad; else; hesh = grad'grad; endif; else; if DOWT; grad = WEIGHT .* grad; endif; hesh = -gradp( &gpnr, b ); @ WEIGHT done in &gpnr @ endif; sgrad = sumc(grad); @ Select only the variables that we want to maximize over @ if CZVAR; sgrad = selif( sgrad, bactive ); hesh = selif( hesh, bactive ); hesh = selif( hesh', bactive ); endif; if (det(hesh)==0); print "Singular Hessian!"; print "Program terminated."; stop; else; ihesh = inv(hesh); direc = ihesh * sgrad; endif; _tol = direc'sgrad; if CZVAR; direc = expand( direc, bactive); endif; print " Tolerance: " _tol; print "--------------------------------------------------------------------------"; if PRTHESS; if CEVAR and CZVAR; printer = expand(sgrad./NP,bactive)~expand(diag(hesh),bactive); printer = nb~EQMAT'*b~EQMAT'*printer[.,1]~EQMAT'printer[.,2]; elseif CEVAR; printer = nb~EQMAT'*b~(EQMAT'*sgrad./NP)~(EQMAT'*diag(hesh)); elseif CZVAR; printer = nb~b~expand(sgrad./NP,bactive)~expand(diag(hesh),bactive); else; printer = nb~b~(sgrad./NP)~(diag(hesh)); endif; print " Coefficients Rel. Grad. Hessian"; else; if CEVAR and CZVAR; printer = expand(sgrad./NP,bactive); printer = nb~EQMAT'*b~EQMAT'*printer[.,1]; elseif CEVAR; printer = nb~EQMAT'*b~(EQMAT'*sgrad./NP); elseif CZVAR; printer = nb~b~expand(sgrad./NP,bactive); else; printer = nb~b~(sgrad./NP); endif; print " Coefficients Rel. Grad."; endif; print printer; if (_tol >= 0) and (_tol < 1.e-6); break; elseif _tol < 0; direc = -direc; endif; step1 = STEP; lambda = .5; repeat = 1; first = 1; _f1 = _f0; steploop: step1 = step1 * lambda; _fold = _f1; if DOWT; _f1 = sumc(WEIGHT .* f(b+step1*direc,XMAT)); else; _f1 = sumc( f( b+step1*direc, XMAT ) ); endif; print "--------------------------------------------------------------------------"; print " Step: " step1; print " Function: " _f1; if repeat; print " Change: " _f1-_f0; else; print " Change: " _f1-_fold; endif; if MNR; if (step1 < 1.e-5); print "Failed to find increase."; retp(b); elseif (_f1 <= _f0) and (repeat); first = 0; goto steploop; elseif (_f1 > _fold) and (first); lambda = 2; repeat = 0; goto steploop; endif; else; if (step1 < 1.e-5); print "Failed to find increase."; retp(b); elseif (_f1 <= _f0); goto steploop; endif; endif; if (repeat); b = b + step1*direc; _f0 = _f1; else; b = b + .5 * step1 * direc; _f0 = _fold; endif; endo; print "=========================================================================="; print; format /m1 /rdt 1,8; if (_tol< EPS); print "Converged with tolerance: " _tol; print "Function value: " _f0; else; print "Stopped with tolerance: " _tol; print "Function value: " _f1; endif; print; lambda = eig(hesh); if lambda>0; print "Function is concave at stopping point."; else; print "WARNING: Function is not concave at stopping point."; endif; print; if ROBUST == 0 and DOWT == 1; @ Create Hessian and gradient for regular @ @ standard errors when have weights. @ if CZVAR; grad = grad'grad; grad = selif(grad, BACTIVE); grad = selif(grad', BACTIVE); ihesh = ihesh*(grad)*ihesh; else; ihesh = ihesh*(grad'grad)*ihesh; endif; endif; if ROBUST == 1; @Create gradient and hessian for robust standard errors. @ if METHOD == 1; hesh = - gradp( &gpnr, b ); if CZVAR; grad = selif( grad', bactive )'; hesh = selif( hesh, bactive ); hesh = selif( hesh', bactive ); endif; endif; ihesh=inv(hesh); ihesh=ihesh*(grad'grad)*ihesh; print "Uses robust standard errors."; else; print "Uses non-robust standard errors."; endif; if CEVAR and CZVAR; printer = expand(sqrt(diag(ihesh)), BACTIVE); printer = nb~EQMAT'*b~EQMAT'*printer; elseif CEVAR; printer = nb~EQMAT'*b~(EQMAT'*sqrt(diag(ihesh))); elseif CZVAR; printer = nb~b~expand(sqrt(diag(ihesh)),BACTIVE); else; printer = nb~b~sqrt((diag(ihesh))); endif; format /m1 /rdt 16,8; print " Parameters Estimates Standard Errors"; print "--------------------------------------------------------------------------"; print printer; retp( b ); endp; /* This proc checks the inputs if VERBOSE=1 */ proc (0) = pcheck; local i,j; @ Checking XMAT @ if (rows(XMAT) /= NOBS); print "XMAT has" rows(XMAT) "rows"; print "but it should have NOBS=" NOBS "rows."; print "Program terminated."; stop; elseif(cols(XMAT) /= (NVAR*NALT)); print "XMAT has" cols(XMAT) "columns"; print "but it should have NVAR*NALT= " (NVAR*NALT) "columns."; print "Program terminated."; stop; else; print "XMAT has:"; print "Rows: " NOBS; print "Cols: " NVAR*NALT ; print "Containing" NVAR "variables for " NALT " alternatives."; print; endif; @ Check dependent variable @ if (rows(YVEC) /= NOBS); print "YVEC has" rows(YVEC) "rows"; print "but it should have NOBS=" NOBS "rows."; print "Program terminated."; stop; endif; if (cols(YVEC) /= 1); print "YVEC should have only one column."; print "Program terminated."; stop; endif; @ Check TIMES @ if (rows(TIMES) /= NP); print "TIMES has" rows(TIMES) "rows"; print "but it should have NP=" NP "rows."; print "Program terminated."; stop; endif; if (cols(TIMES) /= 1); print "TIMES should have only one column."; print "Program terminated."; stop; endif; if (sumc(TIMES) /= NOBS); print "sumc of TIMES =" sumc(TIMES) "but it should equal NOBS."; print "Program terminated."; stop; endif; @ Check WEIGHT @ if (DOWT /=0) and (DOWT /=1); print "DOWT must be 0 or 1."; print "Program terminated."; stop; endif; if DOWT; if (rows(WEIGHT) /= NP); print "WEIGHT has" rows(WEIGHT) "rows"; print "but it should have NP=" NP "rows."; print "Program terminated."; stop; endif; if (cols(WEIGHT) /= 1); print "WEIGHT should have only one column."; print "Program terminated."; stop; endif; endif; @ Checking fixed coefficients @ if (NFC /= 0); if (rows(IDFC) /= NFC); print "IDFC has" rows(IDFC) "rows when it should have NFC=" NFC "rows."; print "Program terminated."; stop; endif; if (cols(IDFC) /= 1); print "Commas are needed between the elements of IDFC."; print "Program terminated."; stop; endif; if (sumc(IDFC.>NVAR)>0); print "IDFC identifies a variable that is not in the data set."; print "All elements of IDFC should be <= NVAR, which is " NVAR; print "Program terminated."; stop; endif; endif; @ Checking normal coefficients @ if (NNC /= 0); if (rows(IDNC) /= NNC); print "IDNC has" rows(IDNC) "rows when it should have NNC=" NNC "rows."; print "Program terminated."; stop; endif; if (cols(IDNC) /= 1); print "Commas are needed between the elements of IDNC."; print "Program terminated."; stop; endif; if (sumc(IDNC.>NVAR)>0); print "IDNC identifies a variable that is not in the data set."; print "All elements of IDNC should be <= NVAR, which is " NVAR; print "Program terminated."; stop; endif; endif; @ Check uniformly distributed coefficients @ if (NUC /= 0); if (rows(IDUC) /= NUC); print "IDUC has" rows(IDUC) "rows when it should have NUC=" NUC "rows."; print "Program terminated."; stop; endif; if (cols(IDUC) /= 1); print "Commas are needed between the elements of IDUC."; print "Program terminated."; stop; endif; if (1-(IDUC <= NVAR)); print "IDUC identifies a variable that is not in the data set."; print "All elements of IDUC should be <= NVAR, which is " NVAR; print "Program terminated."; stop; endif; endif; @ Check traingularly distributed coefficients @ if (NTC /= 0); if (rows(IDTC) /= NTC); print "IDTC has" rows(IDTC) "rows when it should have NTC=" NTC "rows."; print "Program terminated."; stop; endif; if (cols(IDTC) /= 1); print "Commas are needed between the elements of IDTC."; print "Program terminated."; stop; endif; if (1-(IDTC <= NVAR)); print "IDTC identifies a variable that is not in the data set."; print "All elements of IDTC should be <= NVAR, which is " NVAR; print "Program terminated."; stop; endif; endif; @ Checking log-normal coefficients @ if (NLC /= 0); if (rows(IDLC) /= NLC); print "IDLC has" rows(IDLC) "rows when it should have NLC=" NLC "rows."; print "Program terminated."; stop; endif; if (cols(IDLC) /= 1); print "Commas are needed between the elements of IDLC."; print "Program terminated."; stop; endif; if (sumc(IDLC.>NVAR)>0); print "IDLC identifies a variable that is not in the data set."; print "All elements of IDLC should be <= NVAR, which is " NVAR; print "Program terminated."; stop; endif; endif; @ Check CENSOR @ if ((CENSOR /= 1) and (CENSOR /= 0)); print "CENSOR must be 0 or 1."; print "Program terminated."; stop; endif; if CENSOR; if IDCENSOR > NVAR; print "The censoring variable cannot be located"; print "since you set CENSOR larger than NVAR."; print "Program terminated."; stop; endif; i = (IDCENSOR-1)*NALT; j = 1; do while j<=NOBS; if (XMAT[j,i+YVEC[j,1]]==0); print "Your censoring variable eliminates the chosen alternative"; print "for observation " j; print "Program terminated."; stop; endif; j = j + 1; endo; endif; @ Check PRTHESS and RESCALE @ if ((PRTHESS /= 1) and (PRTHESS /= 0)); print "PRTHESS must be 0 or 1."; print "Program terminated."; stop; endif; if ((RESCALE /= 1) and (RESCALE /= 0)); print "RESCALE must be 0 or 1."; print "Program terminated."; stop; endif; @ Check B @ if (rows(B) /= NEVAR); print "Starting values B has " rows(B) "rows"; print "when it should have NFC+2*NNC+2*NUC+2*NTC+2*NLC= " NEVAR "rows."; print "Program terminated."; stop; endif; if (cols(B) /= 1); print "Commas needed between the elements of B."; print "Program terminated."; stop; endif; @ Check CZVAR @ if ((CZVAR /= 1) and (CZVAR /= 0)); print "CZVAR must be 0 or 1."; print "Program terminated."; stop; endif; if CZVAR; if (rows(BACTIVE) /= NEVAR); print "BACTIVE has " rows(BACTIVE) "rows"; print "when it should have NFC+2*NNC+2*NUC+2*NTC+2*NLC= " NEVAR "rows."; print "Program terminated."; stop; endif; if (cols(BACTIVE) /= 1); print "Commas needed between the elements of BACTIVE."; print "Program terminated."; stop; endif; endif; @ Check CEVAR @ if ((CEVAR /= 1) and (CEVAR /= 0)); print "CEVAR must be 0 or 1."; print "Program terminated."; stop; endif; if CEVAR; if (cols(EQMAT) /= NEVAR); print "EQMAT has " cols(EQMAT) " columns"; print "when it should have NFC+2*NNC+2*NUC+2*NTC+2*NLC=" NEVAR " columns."; print "Program terminated."; stop; endif; if (rows(EQMAT)>=NEVAR); print "EQMAT has " rows(EQMAT) " rows"; print "when it should have strictly less than NFC+2*NNC+2*NUC+2*NTC+2*NLC=" NEVAR; print "rows."; print "Program terminated."; stop; endif; endif; @ Checking NREP @ if (NREP <= 0); print "Error in NREP: must be positive."; print "Program terminated."; stop; endif; @ Check SEED1 @ if (SEED1<=0); print "SEED1 =" SEED1 "must be a positive integer."; print "Program terminated."; stop; endif; @ Check METHOD and OPTIM @ if (METHOD < 1); print "METHOD must be 1-6."; print "Program terminated."; stop; endif; if (OPTIM /= 1) and (OPTIM /= 2); print "OPTIM must be 1 or 2."; print "Program terminated."; stop; endif; if ((OPTIM == 1) and (METHOD > 2)); print "Method " METHOD " is not an option with OPTIM = 1."; print "Program terminted."; stop; endif; if ((OPTIM == 2) and (METHOD > 6)); print "Method " METHOD " is not an option with OPTIM = 2."; print "Program terminated."; stop; endif; @ Check MNR @ if ((MNR /= 1) and (MNR /= 0)); print "MNR must be 0 or 1."; print "Program terminated."; stop; endif; @ Check STEP @ if (STEP<=0); print "STEP must be greater than zero."; print "Program terminated."; stop; endif; @ Check ROBUST @ if ((ROBUST /= 1) and (ROBUST /= 0)); print "ROBUST must be 0 or 1."; print "Program terminated."; stop; endif; @ Check DRAWS @ if ((DRAWS /= 1) and (DRAWS /= 2)); print "DRAWS must be 1 or 2."; print "Program terminated."; stop; endif; @ Check HALT @ if ((HALT /= 1) and (HALT /= 2)); print "HALT must be 1 or 2."; print "Program terminated."; stop; endif; print "The inputs pass all the checks and seem fine."; print; retp; endp; /* Halton procedure */ /* Proc to create Halton sequences using the pattern described in Train, "Halton Sequences for Mixed Logit." The integer n is the length of the Halton sequence that is required, and s is the prime that is used in creating the sequence. Given the length of the sequence n, the integer k is determined such that s^k 0.5 and p >= 0.5 @ maskgt = (p .> 0.5); maskeq = (p .ne 0.5); sgn = missrv(miss(maskgt,0),-1); @ Convert p > 0.5 to 1-p @ pn = (maskgt-p).*sgn+mask1+mask0; clear maskgt; @ Computation of function for p < 0.5 @ y=sqrt(abs((-2*ln(pn)))); clear pn; norms = y + ((((y*p4+p3).*y+p2).*y+p1).*y+p0)./ ((((y*q4+q3).*y+q2).*y+q1).*y+q0); clear y; @ Convert results for p > 0.5 and p = 0.5 @ norms=((norms.*sgn).*maskeq).*(1-mask0).*(1-mask1)+mask0.*inf0+mask1.*inf1; retp(norms); endp;