nlog_rd.ado

是一个经济学管理应用软件 很难找的 但是经济学学生又必须用到

*! version 7.0.0  12jul2000

program define nlog_rd 
	version 7.0
	local i 1
	while `i' <= ($nlog_l + $nlog_t) {
		local clist `clist' c`i'
		local i = `i' + 1
	}
	args todo b lnf g negH `clist' 
						/* main eqs */
	local i 1
	while `i' <= $nlog_l {
		tempvar beta`i' 
		mleval `beta`i'' = `b', eq(`i')
		local i = `i' + 1
	}
						/* taus */
	local i 1
	local k = $nlog_l+1 
	while `i' <= $nlog_l - 1 {
		tempvar tau`i'
		qui gen double `tau`i'' = 0 
		local j 1
		while `j' <= $nlog_T[1,`i'] {
			tempvar tau`i'`j' 
			mleval `tau`i'`j'' = `b', eq(`k') 
			local n = (`i'+1)*2
			qui replace `tau`i'' = `tau`i'`j'' /*
			*/ if  ${ML_y`n'} == `j'  
			local k = `k' + 1
			local j = `j' + 1
		}
		local i = `i' + 1
	}
	local bylist$nlog_l $nlog_id 
	local i =$nlog_l -1
	while `i' > 0 {
		local j = `i' + 1
		local n = `j' * 2
		local bylist`i' `bylist`j'' ${ML_y`n'}
		local i = `i' - 1
	}

						/* level 1 */
	tempvar I1 p1
	qui bysort `bylist1': gen double `I1' = sum(exp(`beta1'))
	qui by `bylist1': replace `I1' = `I1'[_N]
	qui gen double `p1' = exp(`beta1')/`I1'
	qui replace `I1' = ln(`I1')
  						/* rest levels */
	local i 2 
	while `i' <= $nlog_l {
		tempvar I`i' p`i'
		local j = `i' - 1
		qui gen double `p`i'' = exp(`beta`i'' + `tau`j''*`I`j'')
		tempvar tmp
		qui gen double `tmp' = 0
		qui bysort `bylist`j'' : replace `tmp' = `p`i'' if _n == 1 
		qui by `bylist`i'': gen double `I`i'' = sum(`tmp')	
		qui by `bylist`i'' : replace `I`i'' = `I`i''[_N]
		qui replace `p`i'' = `p`i''/`I`i''	
		qui replace `I`i'' = ln(`I`i'')
		local i = `i' + 1
	}
						/* ln(p) = sum(ln(p`i')) */
	local i 1
	tempvar lnp
	qui gen double `lnp' = 0
	while `i' <= $nlog_l {
		qui replace `lnp' = `lnp' + ln(`p`i'')  
		local i = `i' + 1
	}	
	qui replace `lnp' = 0 if $ML_y1 == 0
	
	mlsum `lnf' =  `lnp'
	if `todo' == 0 | `lnf' == . { exit }

	if $nlog_l > 3 { exit }
		
						/* BEGIN getting gradients */
						/* main eqs */
	local i 1
	while `i' <= $nlog_l {
		tempvar d`i' temp touse
		qui gen `touse' = 1
		if `i' > 1 {
			local k = `i' - 1
			qui bysort `bylist`k'': replace `touse' = 0 if _n != 1
		}
		tempname g`i'
		qui gen double `temp' = 0
		qui gen double `d`i'' = 0
		local j = `i'
		while `j' <= $nlog_l {
			local m = (`j'-1)*2 + 1
			if `j' == $nlog_l {
				qui replace `temp' =/*
				*/ (${ML_y`m'} - `p`j'')*`touse'
			}
			else {
				local k = 2*`j' + 1
				qui replace `temp' = /*
				*/ (${ML_y`m'}- `p`j''*${ML_y`k'})*`touse'
			}
			local k = `j' - `i'  
			while `k' >= 1 {
				local m = `k' + `i' - 1	
				qui replace `temp' = `temp'* `p`m''*`tau`m''
				local k=`k' - 1
			}
			qui replace `d`i'' = `d`i'' + `temp'	
			local j = `j' + 1
		}
		qui replace `c`i'' = `d`i''
		mlvecsum `lnf' `g`i'' = `d`i'', eq(`i')
		local i = `i' + 1
	}
							/* taus */

	local k = $nlog_l 
	if $nlog_l == 2 {				/* 2 levels */
		local j 1
		while `j' <= $nlog_T[1,1] {
			local k = `k' + 1
			tempname g`k'
			tempvar d1`j' I1`j' p1`j'
			qui by $nlog_id : gen double `I1`j'' = /*
			*/ `I1' if $ML_y4 == `j' 
			qui bysort $nlog_id (`I1`j'') :        /*
			*/ replace `I1`j'' = `I1`j''[_n-1] /*
			*/ if `I1`j'' == . 
			qui by $nlog_id : gen double `p1`j'' = /*
			*/ `p2' if $ML_y4 == `j'
			qui bysort $nlog_id (`p1`j'') :        /*
			*/ replace `p1`j'' = `p1`j''[_n-1]  if `p1`j'' == . 
			qui gen double `d1`j'' = $ML_y1 *   /*
			*/ `I1`j''*(cond($ML_y4==`j', 1, 0)- `p1`j'' )
			qui replace `c`k'' = `d1`j''
			mlsum `g`k'' = `d1`j''
			local j = `j' + 1
		}
	}
	if $nlog_l == 3 {				/* 3 levels */
		  local j 1 
		  while `j' <= $nlog_T[1,1] {
			local k = `k' + 1
			tempname g`k'
			tempvar d1`j' d2`j' I1`j' I2`j' p1`j' p2`j' touse
							
						/* dlnp2/dtau1  */
			qui gen `touse' = $ML_y6 if $ML_y4 == `j' 
			qui sort `touse'
			local t1 = `touse'[1]
			qui bysort `bylist2' : gen double `I1`j'' = /*
			*/ `I1' if $ML_y4 == `j' 
			qui bysort $nlog_id (`I1`j'') :             /*
			*/ replace `I1`j'' = `I1`j''[_n-1]  if `I1`j'' == . 
			qui bysort `bylist2' : gen double `p1`j'' = /*
			*/ `p2' if $ML_y4 == `j'
			qui bysort $nlog_id (`p1`j''):              /*
			*/ replace `p1`j'' = `p1`j''[_n-1] if `p1`j'' == . 
			qui gen double `d1`j'' =cond($ML_y6==`t1',1,0)*    /*
			*/ $ML_y1 * `I1`j''*(cond($ML_y4==`j', 1, 0)- `p1`j'' )
			
						/* dlnp3/dtau1 */	
			qui by $nlog_id : gen double `p2`j'' = /*
			*/ `p3' if $ML_y6 == `t1'
			qui bysort $nlog_id (`p2`j'') :             /*
			*/ replace `p2`j'' = `p2`j''[_n-1]  if `p2`j'' == . 
			qui gen double `d2`j'' = `tau2`t1''*$ML_y1 *   /*
			*/ `I1`j''*(cond($ML_y6==`t1', 1, 0)- `p2`j'')*`p1`j''
		
			qui replace `c`k'' = `d1`j'' + `d2`j''	
			mlsum `g`k'' = `c`k'' 
			local j = `j' + 1
		}

						/* dlnp3/dtau2 */
		local j = 1
		while `j' <= $nlog_T[1,2] {     
			local k = `k' + 1
			tempname g`k'
			tempvar d2`j' I2`j' p2`j'
			qui by $nlog_id : gen double `I2`j'' = /*
			*/ `I2' if $ML_y6 == `j' 
			qui bysort $nlog_id (`I2`j''):         /*
			*/ replace `I2`j'' = `I2`j''[_n-1] if `I2`j'' == . 
			qui by $nlog_id : gen double `p2`j'' = /*
			*/ `p3' if $ML_y6 == `j'
			qui bysort $nlog_id (`p2`j'') :        /*
			*/  replace `p2`j'' = `p2`j''[_n-1] if `p2`j'' == . 
			qui gen double `d2`j'' = $ML_y1 *   /*
			*/ `I2`j''*(cond($ML_y6==`j', 1, 0)- `p2`j'' )
			qui replace `c`k'' = `d2`j''
			mlsum `g`k'' = `c`k''
			local j = `j' + 1
		}
	}
			 	
			 
	mat `g' = `g1'
	local i 2
	while `i' <= ($nlog_t+$nlog_l) {
		mat `g' = (`g', `g`i'')
		local i = `i' + 1
	}	
end	

exit

Example 3 levels : (refer to Limdep manual page 591-592)

		math				code

p        = p(k|i,j)*p(j|i)*p(i)			p  = p1*p2*p3

I(ij)    = sum(p(k|i,j))			I1 = sum(p1), by(2,3)
	    k
I(i)     = sum(p(j|i))				I2 = sum(p2), by(3)	
       	    j	

p(k|i,j) = exp(beta1*x)/exp(I(ij))               p1 = exp(`beta1')/exp(I1)

p(i|j)   = exp(beta2*y + tauij*I(ij))/exp(I(i))  p2 = exp(`beta2' + tau1*I1)/..

p(i)     = exp(beta3*z + taui*I(i))/sum(...)     p3 = exp(`beta3' + tau2*I2)/.. 


gradiants (use mlvecsum):

level 1:  use the whole dataset

dlnp1/db1 = (ML_y1 - p1*ML_y3)        	   
dlnp2/db1 = (ML_y3 - p2*ML_y5)*tau1*p1     
dlnp3/db1 = (ML_y5 - p3)*tau1*p1*tau2*p2

level 2: keep one obs for each level of ML_y3 within each id

dlnp2/db2 = (ML_y3 - p2*ML_y5)
dlnp3/db2 = (ML_y5 - p3)*tau2*p2

level 3: keep one obs for each level of ML_y5 within each id 

dlnp3/db3 = (ML_y5 - p3)

taus:

dlnp(j|i)/dtau(m|l) = 1(l=i)[1(m=j) - p(m|l)]*I(m|l)

dlnp(i)/dtau(m|l) = tau(i)[1(l=i) - p(l)]*p(m|l)*I(m|l)

dlnp(i)/dtau(l) = [1(l=i) - p(l)]*I(l)