/* Bas van der Klaauw */ /* Date: 10-28-1998 */ /* This GAUSS-program estimates the sample selection model using */ /* a flexble specification for the joint distribution of the error */ /* terms. This flexible specification is based on Hermite series */ /* with the bivariate normal density as the base function, which */ /* straightforwardly allows for testing the normality assumption. */ /* This flexible parametric density function and the normality test */ /* are discussed in Bas van der Klaauw and Ruud H. Koning (2000), */ /* `Testing the normality assumption in the Sample Selection Model */ /* with an Application to Travel Demand. The flexible parametric */ /* density is motived on Gallant and Nychka (1987), `Semi-nonparametric */ /* Maximum Likelihood Estimation', Econometrica 55, p 363-390. */ new; n1 = ..; /* Number of exogenous variables in the regression equation */ n2 = ..; /* Number of exogenous variables in the choice equation */ nn = ..; /* Number of observations */ kk = ..; /* K */ load data[nn,(1+n1+n2)] = ..; /* Data file should consist the endogenous variable, the exogenous */ /* variables in the regression equation and the exogenous variables */ /* in the choice equation. */ library cml; #include cml.ext; cmlset; _cml_EqProc = ℰ _cml_GradTol = 1e-4; _cml_Algorithm = 1; _cml_CovPar = 2; _cml_MaxIters = 200; output file = samgen.out reset; startv = ......; /* Startvalues of the optimization, this consist of startvalues for */ /* the parameters in the regression equation, parameters in the choice */ /* equation, delta, sigma and the alpha's. */ mom = momp; {schatpar, f, g, cov, retcode} = cmlprt(cml(data,0,&loglik,startv)); output off; /* Loglikelihood function */ /* Assumed are the identifying restrictions */ /* alpha[0,0] = 1 (normalization of the alpha's) */ /* The expectations of the error terms are fixed (location of the density function) */ /* delta2 = sqrt(2) (scale of the choice equation) */ proc loglik(param, data); local s, i, j, k, l, hulp, logl, beta1, beta2, alpha, delta, sigma, w1, w2, y, x1, x2, z, p1, p2, p3, b1, b2; b1 = param[1]; beta1 = param[2:n1+1]; b2 = param[n1+2]; beta2 = param[n1+3:n1+n2+2]; delta = param[n1+n2+3]; sigma = param[n1+n2+4]; alpha = zeros(kk+1,kk+1); i = 0; do while i <= kk; j = 0; do while j <= kk; alpha[i+1,j+1] = param[n1+n2+4+j+i*(kk+1)]; j = j + 1; endo; i = i + 1; endo; alpha[1,1] = 1; y = data[.,1]; x1 = data[.,2:n1+1]; x2 = data[.,n1+2:n1+n2+1]; w1 = y-x1*beta1-b1; w2 = -x2*beta2-b2; z = (y.==0); p1 = int1(w1,w2,delta,sigma); p2 = int2(w2,delta,sigma); p3 = int3(delta,sigma); hulp = zeros(nn,1); s = 0; i = 0; do while i <= kk; j = 0; do while j <= kk; k = 0; do while k <= kk; l = 0; do while l <= kk; s = s + alpha[i+1,j+1]*alpha[k+1,l+1]*p3[i+k+1,j+l+1]; hulp = hulp + (1-z).*(alpha[i+1,j+1]*alpha[k+1,l+1]* p1[.,(2*kk+1)*(i+k)+(j+l)+1]); hulp = hulp + z.*(alpha[i+1,j+1]*alpha[k+1,l+1]* p2[.,(2*kk+1)*(i+k)+(j+l)+1]); l = l + 1; endo; k = k + 1; endo; j = j + 1; endo; i = i + 1; endo; logl = ln(hulp) - ln(s); retp(logl); endp; /* Other functions */ proc expectation(param); local s, i, j, k, l, hulp, logl, beta1, beta2, alpha, delta, sigma, w1, w2, y, x1, x2, z, p1, p2, p3, b1, b2, v1, v2, lambda1, lambda2; delta = param[n1+n2+3]; sigma = param[n1+n2+4]; alpha = zeros(kk+1,kk+1); i = 0; do while i <= kk; j = 0; do while j <= kk; alpha[i+1,j+1] = param[n1+n2+4+j+i*(kk+1)]; j = j + 1; endo; i = i + 1; endo; alpha[1,1] = 1; p3 = int3(delta,sigma); hulp = zeros(nn,1); s = 0; v1 = 0; v2 = 0; i = 0; do while i <= kk; j = 0; do while j <= kk; k = 0; do while k <= kk; l = 0; do while l <= kk; s = s + alpha[i+1,j+1]*alpha[k+1,l+1]*p3[i+k+1,j+l+1]; v1 = v1 + alpha[i+1,j+1]*alpha[k+1,l+1]*p3[i+k+2,j+l+1]; v2 = v2 + alpha[i+1,j+1]*alpha[k+1,l+1]*p3[i+k+1,j+l+2]; l = l + 1; endo; k = k + 1; endo; j = j + 1; endo; i = i + 1; endo; retp((v1/2)|(v2/s)); endp; proc ind(k,delta); local hulp; if k; if k-1; hulp = (k-1)/2*delta^2*ind(k-2,delta); else; hulp = 0; endif; else; hulp = delta * sqrt(pi); endif; retp(hulp); endp; proc kans1(k,onder); local hulp; if k; if k-1; hulp = onder.^(k-1).*exp(-0.5*onder.^2)+(k-1)*kans1(k-2,onder); else; hulp = exp(-0.5*onder.^2); endif; else; hulp = sqrt(2*pi)*cdfnc(onder); endif; retp(hulp); endp; proc kans2(k,boven); local hulp; if k; if k-1; hulp = -boven.^(k-1).*exp(-0.5*boven.^2)+(k-1)*kans2(k-2,boven); else; hulp = -exp(-0.5*boven.^2); endif; else; hulp = sqrt(2*pi)*cdfn(boven); endif; retp(hulp); endp; proc momp; local i, j, k, mm; k = 2*kk+2; mm = zeros(k,k); i = 1; mm[1,1] = 1; do while i <= k-1; mm[i+1,1] = 1; j = 1; do while j <= i/2; mm[i+1,j+1] = mm[i,j+1] + mm[i,j] * (i-1-2*(j-1)); j = j + 1; endo; i = i + 1; endo; retp(mm); endp; proc int3(delta,sigma); local integral, i, j, k, hulp, mu, var, rho, x; mu = sigma; rho = sigma/delta; var = sqrt(delta^2*(1-rho^2)); integral = zeros(2*kk+2,2*kk+2); x = zeros(4*kk+3,1); i = 0; do while i <= 4*kk+2; x[i+1] = ind(i,sqrt(2)); i = i + 1; endo; i = 0; do while i <= 2*kk+1; j = 0; do while j <= 2*kk+1; hulp = 0; k = 0; do while k <= i; hulp = hulp + mom[i+1,k/2+1] * mu^(i-k) * var^k * x[j+i-k+1] / sqrt(2*pi); k = k + 2; endo; integral[i+1,j+1] = hulp; j = j + 1; endo; i = i + 1; endo; retp(integral); endp; proc int2(b,delta,sigma); local integral, i, j, k, hulp, mu, var, rho, x; mu = sigma; rho = sigma/delta; var = sqrt(delta^2*(1-rho^2)); integral = zeros(nn,(2*kk+1)^2+1); x = zeros(nn,4*kk+1); i = 0; do while i <= 4*kk; x[.,i+1] = kans2(i,b); i = i + 1; endo; i = 0; do while i <= 2*kk; j = 0; do while j <= 2*kk; hulp = zeros(nn,1); k = 0; do while k <= i; hulp = hulp + mom[i+1,k/2+1] * mu^(i-k) * var^k * x[.,j+i-k+1] / sqrt(2*pi); k = k + 2; endo; integral[.,(2*kk+1)*i+j+1] = hulp; j = j + 1; endo; i = i + 1; endo; retp(integral); endp; proc int1(a,b,delta,sigma); local integral, i, j, k, hulp, mu, var, rho, x, onder; mu = sigma/delta^2*a; rho = sigma/delta; var = sqrt((1-rho^2)); onder = (b-mu)/var; integral = zeros(nn,(2*kk+1)^2+1); x = zeros(nn,4*kk+1); i = 0; do while i <= 4*kk; x[.,i+1] = kans1(i,onder); i = i + 1; endo; i = 0; do while i <= 2*kk; j = 0; do while j <= 2*kk; hulp = zeros(nn,1); k = 0; do while k <= j; hulp = hulp + (j!/(k!*(j-k)!) * (mu/var).^(j-k)) .* (x[.,k+1] / sqrt(2*pi)); k = k + 1; endo; integral[.,(2*kk+1)*i+j+1] = a^i .* (1/delta * pdfn(a/delta)) .* (var^j * hulp); j = j + 1; endo; i = i + 1; endo; retp(integral); endp;