_vecu.ado

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

			t(name)			///
			s00(name)		///
			s01(name)		///
			s11(name)		///
			val(name)		///
			k(integer)
	
	tempname pi mat1 omega2 temp omegai

	mat `pi' = `alpha'*`beta''
	mat `mat1' = `pi'*(`s01'')
	
	GetOmega , 	s00(`s00')		///
			s01(`s01')		///
			s11(`s11')		///
			alpha(`alpha')		///
			beta(`beta')		///
			omega(`omega2')		

	scalar `val'  = -.5*`t'*`k'*ln(2*_pi) -		///
		.5*`t'*ln(abs(det(`omega'))) -.5*`t'*`k'

end

program define GetAlpha

	syntax , 	gmat(name)		///
			omegai(name)		///
			beta(name)		///
			s11(name)		///
			vecpils(name)		///
			alpha(name)

	tempname psi vec_alphap
	if "`gmat'" == "NONE"  {
		mat `psi' = syminv(`omegai'#(`beta''*`s11'*`beta'))*	///
			(`omegai'#(`beta''*`s11'))*`vecpils'
// in this case psi = vec_alphap
		va2a , vec_alphap(`psi') alpha(`alpha')
	}
	else {
		mat `psi' = syminv(`gmat''*(`omegai'#(`beta''*`s11'*`beta')) ///
			*`gmat')*`gmat''*(`omegai'#(`beta''*`s11'))*`vecpils'
		mat `vec_alphap' = `gmat'*`psi'	
		va2a , vec_alphap(`vec_alphap') alpha(`alpha')

	}


end

program define GetOmega 
	syntax , 	s00(name)		///
			s01(name)		///
			s11(name)		///
			alpha(name)		///
			beta(name)		///
			omega(name)		

	tempname mat1 pi

	mat `pi' = `alpha'*`beta''

	mat `mat1' = `s01'*(`pi'')

	mat `omega' = `s00' - `mat1' - `mat1'' + `pi'*`s11'*`pi''

end

program define GetBeta

	syntax , 	hmat(name)		///
			alpha(name)		///
			omegai(name)		///
			s11(name)		///
			vecpils(name)		///
			ip1(name)		///
			beta(name)		///
			[			///
			amat(name)		///
			]

	tempname phi vec_beta
	if "`hmat'" == "NONE" {
		mat `phi' = syminv((`alpha''*`omegai'*`alpha')#`s11')* 	///
			((`alpha''*`omegai')#`s11')*`vecpils'
// in this case phi = vec(beta)

		vb2b , vecbeta(`phi') beta(`beta')

	}
	else {

		mat `phi' = syminv(`hmat''*((`alpha''*`omegai'*`alpha')#  ///
			`s11')*`hmat')*`hmat''*((`alpha''*`omegai')#`s11') ///
			*(`vecpils' - (`alpha'#`ip1')*`amat')
		mat `vec_beta'	= `hmat'*`phi' + `amat'
		vb2b , vecbeta(`vec_beta') beta(`beta')
	}
end

program define vb2b

	syntax  , vecbeta(name) 		///
		  beta(name)

	local rows = rowsof(`beta')
	local cols = colsof(`beta')

	tempname beta2

	if rowsof(`vecbeta') != `rows'*`cols' | colsof(`vecbeta') != 1 {
		di as err "parameter estimates of beta do not have "	///
			"correct dimension"
		exit 498
	}

	forvalues i = 1/`cols' {
		mat `beta2' = (nullmat(`beta2'), 		///
			`vecbeta'[(`i'-1)*`rows' + 1..`i'*`rows',1] )
	}

	mat `beta' = `beta2'
end		

program define va2a 
	syntax ,	vec_alphap(name)	///
			alpha(name)
			

	local r = colsof(`alpha')
	local k = rowsof(`alpha')

	tempname alpha2

	if rowsof(`vec_alphap') != `r'*`k' | colsof(`vec_alphap') != 1 {
		di as err "parameter estimates of alpha do not have "	///
			"the correct dimension"
		exit 498	
	}		

	forvalues i = 1/`k' {

		mat `alpha2' = (nullmat(`alpha2') \ 	///
			(`vec_alphap'[(`i'-1)*`r' + 1..`i'*`r', 1])' )
	}

	mat `alpha' = `alpha2'
end

program define 	Getsmats 

	syntax , 			///
		r0t(varlist ts) 	///
		r1t(varlist ts) 	///
		touse(varname)		///
		s00(name) 		///
		s01(name) 		///
		s10(name) 		///
		s11(name)  		///
		t(name) 		///
		k(integer)

	tempname bigs

	qui mat accum `bigs' = `r0t' `r1t' if `touse', nocons
	scalar `t' = r(N)
	mat `bigs' = (1/`t')*`bigs'

	mat `s00' =`bigs'[1..`k',1..`k']
	mat `s01' =`bigs'[1..`k',`k'+1...]
	mat `s10' = `s01''
	mat `s11' =`bigs'[`k'+1...,`k'+1...]

end		


program define Johansen 

	syntax , alpha(name)	///
		beta(name)	///
		omega(name)	///
		lam(name)	///
		s10(name)	///
		s01(name)	///
		s11(name)	///
		s00(name)



	tempname s11_sqrt bigs2 s00i u 
	
	MATsqrt `s11_sqrt' `s11'


	mat `s00i' = syminv(`s00')
	mat `bigs2' = `s11_sqrt'*`s10'*`s00i'*`s01'*`s11_sqrt'
	
	mat symeigen `u' `lam' = `bigs2'
	mat `beta' = (`s11_sqrt'*`u')

	mat `alpha' = `s01'*`beta'

	mat `omega' = `s00' -`alpha'*`alpha''

end

program define Mkceqj

	syntax newvarlist, beta2(name) 

	local j 0
	foreach ceq of local varlist {
		local ++j
		local ceqj : word `j' of `varlist'
		mat score double `ceqj' = `beta2' , eq(#`j')
		global S_madece madece
	}

end


program define Getb

		syntax varlist(ts numeric) 		///
			[if] , 				///
			pm1(integer) 			///
			exfront(varlist ts)		///
			df(integer)			///
			bmat(name)			///
			vmat(name)			///
			dvlist2(namelist)		///
			zbetan(string)			///
			[				///
			exog(varlist ts)		///
			const(numlist)			///
			from(passthru)			///
			quietly				///
			]

		marksample touse

		if "`quietly'" == "" di 

		if "`const'" != "" {
			`quietly' _cvar `varlist' if `touse', 		///
				lags(1/`pm1')				///
				const(`const') exog(`exog')		///
				noconstant exfront(`exfront')		///
				dfmodel(`df') isure tol(1e-9)

// adjustment on next line properly counts df when there are 
// constraints on alpha
			mat `vmat' = (e(N)/(e(N)-`df'))*e(V)
			
		}
		else {
			if `pm1' > 0 {
				local qslags " lags(1/`pm1') "
			}
			else {
				local qslags " zlags "
			}
			`quietly' _qsur `varlist' if `touse', `qslags'	///
				 exfront(`exfront') exog(`exog')	///
				 noconstant dfmodel(`df')
			mat `vmat' = e(V)
		}		

		mat `bmat' = e(b)

		foreach evar of local dvlist2 {
			foreach zvar of local zbetan {
				local vnames "`vnames' `evar':`zvar' "
			}
		}


		mat colnames `bmat' = `vnames'
		mat colnames `vmat' = `vnames'
		mat rownames `vmat' = `vnames'

end			

//  This program calculates the rank of the full jacobian and the rank of the
//  jacobian with respect to beta only for the values of alpha and beta passed
//  or for the ``consecutive full rank" matrix

program define GetDF, rclass

	syntax , gmat(string) 		///
		hmat(string) 		///
		beta(name) 		///
		alpha(name)		///
		k(integer)		///
		nz1(integer)		///
		rank(integer)		///
		[			///
		rmat(string)		///
		]

	tempname rbeta ralpha
	
	local row_alpha = rowsof(`alpha')
	local col_alpha = colsof(`alpha')

	local row_beta = rowsof(`beta')
	local col_beta = colsof(`beta')

	MKcfrm `rbeta' `row_beta' `col_beta'
	MKcfrm `ralpha' `row_alpha' `col_alpha'

	CKjacob , k(`k') beta(`beta') alpha(`alpha') 		///
		gmat(`gmat') hmat(`hmat') nz1(`nz1') 		///
		rank(`rank')

	local jacob_cols     = r(jacob_cols)	
	local jacob_hat_rank = r(jacob_rank)
	local jacob_hat_lrdf = r(df_lr)

	CKjacob_b , k(`k') beta(`rbeta') hmat(`hmat') 		///
		rmat(`rmat') nz1(`nz1') rank(`rank')	

	local beta_iden      = r(beta_iden)
	local beta_iden_cnt  = r(beta_iden_cnt)

	ret scalar jacob_cols      = `jacob_cols'
	ret scalar jacob_hat_rank  = `jacob_hat_rank'
	ret scalar jacob_hat_lrdf  = `jacob_hat_lrdf'

	ret scalar beta_iden       =  `beta_iden'
	ret scalar beta_iden_cnt   =  `beta_iden_cnt'
end

// This program makes and checks the rank of a jacobian for a 
// given alpha and beta

program define CKjacob , rclass

	syntax ,	k(integer) 		///
			beta(name)		///
			alpha(name) 		///
			gmat(string) 		///
			hmat(string)		///
			nz1(integer)		///
			rank(integer)

	tempname j1 j2 jacob idmat

	if "`gmat'" == "NONE" {
		mat `j1' = I(`k')#`beta'
	}
	else {
		mat `j1' = (I(`k')#`beta')*`gmat'
	}

	if "`hmat'" == "NONE" {
		mat `j2' = `alpha' # I(`nz1')
	}
	else {
		mat `j2' = (`alpha' # I(`nz1'))*`hmat'
	}

	mat `jacob' = `j1', `j2'

	mat `idmat' = syminv(`jacob''*`jacob')

	local c1 = colsof(`idmat')
		
	local df_b 0 
	forvalues i = 1/`c1' {
		if `idmat'[`i',`i'] > 5e-16 local ++df_b
	}	


	local df_lr = ( `k' + `nz1' -`rank')*`rank' -`df_b'

	ret scalar jacob_cols = `c1'
	ret scalar jacob_rank = `df_b'

	ret scalar df_lr = `df_lr'

end			

// This program checks the jacobian for beta only for a given beta
program define CKjacob_b , rclass

	syntax ,	k(integer) 		///
			beta(name)		///
			hmat(string)		///
			nz1(integer)		///
			rank(integer)		///
			[			///
			rmat(string)		///
			]


	if "`hmat'" == "NONE" {
		ret scalar beta_iden_cnt = 0
		ret scalar beta_iden     = 0
	}
	else {

		tempname brmat

		mat `brmat' = `rmat'*(I(`rank') # `beta')
		mat `brmat' = syminv(`brmat''*`brmat')
		local brcols = colsof(`brmat')
		local beta_iden_cnt 0
		forvalues i = 1/`brcols' {
			if `brmat'[`i',`i'] > 1e-15 {
				local ++ beta_iden_cnt
			}
		}

		if `beta_iden_cnt' >= `rank'^2 {
			ret scalar beta_iden = 1
		}
		else {
			ret scalar beta_iden = 0
		}

		ret scalar beta_iden_cnt = `beta_iden_cnt'

	}

end

// This program makes a ``consecutive full rank" matrix

program define MKcfrm 
	
	args mat rows cols

	mat `mat' = J(`rows', `cols', 0)

	local cnt 0

	forvalues j=1/`cols' {
		forvalues i=1/`rows' {
			local ++cnt
			mat `mat'[`i',`j'] = `cnt'
		}
	}
	mat `mat' = .01*`mat'


end

program define vec_rmcoll, rclass

	syntax varlist(ts), 		///
		touse(varname)		///
		p(integer)		///
		[			///
		rmcoll			///
		ownonly			///
		]



	_rmcoll `varlist' if `touse'
	local varlist "`r(varlist)'"

	local pm1 = `p' - 1

	if "`ownonly'" == "" {
		Getz2base `varlist', pm1(`pm1')
		local z2base "`r(varlist)'"
	}

	foreach dvar of local varlist {
		if "`ownonly'" != "" {
			tsunab dvar0 : `dvar'
			Getz2base `dvar0', pm1(`pm1')
			local z2base "`r(varlist)'"

		}

		tsunab ddvar : d.`dvar'
		local clean 0
		while `clean' == 0 {
			local ext_vlist `ddvar' `z2base'
			
			if "`rmcoll'" == "" {
				qui reg `ext_vlist' if `touse'
				if reldif(_b[_cons], 0) < 5e-16 {
					local clean 0
				}
				else {
					local clean 1
				}
			}
			else {
				qui _rmcoll `ext_vlist'
				local rm_vlist "`r(varlist)'"
				local clean : list ext_vlist == rm_vlist
			}


			if `clean' == 0 {
				local p = `p' - 1
				if `p' <1 {
					di as err "cannot fit model "	///
						"due to collinearity"
					exit 498
				}

				local pm1 = `p' - 1
				di as txt "maximum lag reduced to "	///
					"`p' due to collinearity"
				local reduced `reduced' `p'	

				if "`ownonly'" == "" {
					Getz2base `varlist', pm1(`pm1')
					local z2base "`r(varlist)'"
				}
				else {
					tsunab dvar0 : `dvar'
					Getz2base `dvar0', pm1(`pm1')
					local z2base "`r(varlist)'"
		
				}
			}
		}
	}

	ret local varlist "`varlist'"
	ret scalar p = `p'
	ret local reduced "`reduced'"
end

program define Getz2base, rclass

	syntax varlist(ts), pm1(integer) 

	if `pm1' > 0 {
		tsunab z2base : l(1/`pm1').d.(`varlist') 
	}
	else {
		local z2base
	}

	ret local varlist "`z2base'"
end

program define DerrMsg
	di as err "cannot perform estimation using Johansen's method"

end
exit