@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ @ @ @ Problem Set #2 - Econ 673 - Spring 2008 @ @ Question 2 @ @ @ @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ @ @ @ Note: MAXLIK is a special routine and you have to tell GAUSS that your @ @ will be using it @ @ @ @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ library maxlik; #include maxlik.ext; maxset; output file = Prob3q2.out reset; @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ @ @ @ The first step in this program is to generate the data @ @ @ @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ /* specify parameter values */ beta0 = -1.0; betag = 1.2; betap = -.004; betac = 1.0; betae = 1.0; beta = beta0|betag|betap|betac|betae; /* generate explanatory variables */ nobs = 1000; a = 2.5; b = 3.5; mu = 2.9; sig = 0.5; ua = cdfn((a - mu)/sig); ub = cdfn((b - mu)/sig); u = ua + (ub - ua)*rndu(nobs,1); g = mu + sig*cdfni(u); /*gasoline prices*/ p = 500 + 1500*rndu(nobs,1); /*hybrid price premium*/ Dc = (rndu(nobs,1).>0.75); /*college graduate*/ De = (rndu(nobs,1).>0.65); /*membership in env. group*/ X = ones(nobs,1)~g~p~Dc~De; k = cols(X); ystar = X*beta - rndn(nobs,1); y = (ystar.>0); q = 2*y - 1; @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ @ @ @ Generate Summary Statistics @ @ @ @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ z = y~X; meanz = meanc(z)'; stdz = stdc(z)'; maxz = maxc(z)'; minz = minc(z)'; print ; print "summary statistics"; print "**********************************************"; print " mean std min max "; format /m1 /rd 10,4; length = cols(meanz); for i (1,length,1); print meanz[1,i] stdz[1,i] minz[1,i] maxz[1,i]; endfor; print ; @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ @ @ @ Do-loop for Newton Raphson routine @ @ @ @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ toler = 0.000001; convg = 0; startv = zeros(k,1); iter = 1; imax = 20; betaold = startv; do while ((iter<=imax).and(convg == 0)); index = X * betaold; @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ @ @ @ Compute gradient @ @ @ @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ lambda = q .* pdfn(index) ./ cdfn(q .* index); grad = (lambda'*X)'; @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ @ @ @ Compute Hessian @ @ @ @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ hess1 = zeros(k,k); for i (1,nobs,1); hess1 = hess1 - lambda[i,1]*(lambda[i,1] + index[i,1]) *X[i,.]'*X[i,.]; endfor; @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ @ @ @ Compute updated parameter vector and check for convergence @ @ @ @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ delta = inv(hess1)*grad; betanew = betaold - delta; reldel = delta ./ betaold; info = betaold~betanew~delta~reldel; if abs(maxc(reldel)) < toler; convg = 1; endif; betaold = betanew; print "interation: " iter "coefficients are" betanew; iter = iter + 1; endo; if convg == 1; print "convergence achieved after " iter "iterations"; else; print "convergence not achieved after " iter "iterations"; endif; beta = betanew; @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ @ @ @ Generate Asymptotic Standard Errors @ @ @ @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ lambda1 = pdfn(index) ./ cdfn(index); lambda0 = -pdfn(index) ./ cdfn(-index); Ehess = zeros(k,k); for i (1,nobs,1); Ehess = Ehess - lambda1[i,1]*lambda0[i,1]*X[i,.]'*X[i,.]; endfor; varcov = inv(Ehess); stderr = sqrt(diag(varcov)); print ; print "results of estimation"; print "***********************************"; print " estimate std. err."; print beta~stderr; @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ @ @ @ Compute fitted probabilities and marginal effects @ @ @ @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ riter = 1000; XA = {1 2.5 1500 0 0}; XB = {1 3.5 500 1 1}; XC = {1 3.0 1000 1 0}; XA1 = {1 2.5 1500 0 1}; XB1 = {1 3.5 500 1 1}; XC1 = {1 3.0 1000 1 1}; X1 = X - (zeros(nobs,2)~(500*ones(nobs,1))~zeros(nobs,2)); pfitA = zeros(riter,1); pfitB = zeros(riter,1); pfitC = zeros(riter,1); mfitA = zeros(riter,1); mfitB = zeros(riter,1); mfitC = zeros(riter,1); dfitA = zeros(riter,1); dfitB = zeros(riter,1); dfitC = zeros(riter,1); rfit = zeros(riter,1); betar = beta' + rndn(riter,5) * chol(varcov); for i (1,riter,1); pfitA[i,1] = cdfn(XA*betar[i,.]'); pfitB[i,1] = cdfn(XB*betar[i,.]'); pfitC[i,1] = cdfn(XC*betar[i,.]'); mfitA[i,1] = pdfn(XA*betar[i,.]')*betar[riter,2]; mfitB[i,1] = pdfn(XB*betar[i,.]')*betar[riter,2]; mfitC[i,1] = pdfn(XC*betar[i,.]')*betar[riter,2]; dfitA[i,1] = cdfn(XA1*betar[i,.]') - pfitA[i,1]; dfitB[i,1] = cdfn(XB1*betar[i,.]') - pfitB[i,1]; dfitC[i,1] = cdfn(XC1*betar[i,.]') - pfitC[i,1]; rfit[i,1] = meanc(cdfn(X1*betar[i,.]')) - meanc(cdfn(X*betar[i,.]')); endfor; print ; print "Fitted probabilities"; print "***********************************"; print " type estimate std. err. "; print " I " meanc(pfitA) stdc(pfitA); print " II " meanc(pfitB) stdc(pfitB); print " III " meanc(pfitC) stdc(pfitC); print ; print "Marginal Impact of Gas Prices "; print "***********************************"; print " type estimate std. err."; print " I " meanc(mfitA) stdc(mfitA); print " II " meanc(mfitB) stdc(mfitB); print " III " meanc(mfitC) stdc(mfitC); print ; print "Change in Pr(y=1) from joining Env. Group"; print "*****************************************"; print " type estimate std. err. "; print " I " meanc(dfitA) stdc(dfitA); print " II " meanc(dfitB) stdc(dfitB); print " III " meanc(dfitC) stdc(dfitC); conf = 0.95; print ; print "Change in Pr(y=1) from $500 rebate"; print "***************************************************"; print " " 100*conf "percent confidence interval"; print "estimate std. err. lower bound upper bound"; srebate = sortc(rfit,1); lb = srebate[floor((1 - conf)*riter/2),1]; ub = srebate[ceil((1 + conf)*riter/2),1]; print meanc(rfit) stdc(rfit) lb ub; @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ @ @ @ Estimate the model using maxlik @ @ @ @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ b0 = beta; b0 = {1,0,0,0,0}; dta = y~X; _max_GradProc = &gradprbt; {b_ml,logl,g,varcov,ret}=maxprt(maxlik(dta,0,&prbt,b0)); stderr = sqrt(diag(varcov)); print ; print "maxlik results of estimation"; print "***********************************"; print " estimate std. err."; print b_ml~stderr; end; @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ @ @ @ Procedure to compute loglikelihood value for probit model @ @ @ @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ proc prbt(b,dta); local yvar, xvar, qvar, llik, nobs, k; yvar = dta[.,1]; qvar = 2 * yvar - 1; nobs = rows(dta); k = cols(dta); xvar = dta[.,2:k]; llik = ln(cdfn(qvar .* xvar * b)); /* Likelihood function for probit */ retp (llik); endp; @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ @ @ @ Procedure to compute gradient of loglikelihood value for probit model @ @ @ @@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@@ proc gradprbt(b,dta); local yvar, xvar, index, qvar, gradp, lambda, nobs, k; yvar = dta[.,1]; qvar = 2 * yvar - 1; nobs = rows(dta); k = cols(dta); xvar = dta[.,2:k]; index = xvar*b; lambda = qvar .* pdfn(index) ./ cdfn(qvar .* index); gradp = lambda.*xvar; retp (gradp); endp;