**! version 1.0.1 10May2020
*! author Gabriele Rovigatti, Bank of Italy, Rome, Italy. mailto: gabriele.rovigatti@gmail.com  |  gabriele.rovigatti@bancaditalia.it

/***************************************************************************
** Stata program for Markup Estimation - custom PRODEST version
** Programmed by:   Gabriele Rovigatti
**************************************************************************/

capture program drop prodest_m
program define prodest_m, sortpreserve eclass
	version 10.0
	syntax varlist(numeric min=1 max=1) [if] [in], free(varlist numeric min=1) /*
		  */ proxy(varlist numeric min=1 max=2) state(varlist numeric min=1) /*
		  */ [control(varlist min=1) ENDOgenous(varlist min=1) id(varlist min=1 max=1) t(varlist min=1 max=1) reps(integer 5) /*
		  */ VAlueadded Level(int ${S_level}) OPTimizer(namelist min=1 max=3) MAXiter(integer 10000) /* OPTIMIZER: THERE IS THE POSSIBILITY TO WRITE technique(nr 100 nm 1000) with different optimizers after the number of iterations.
		  */ poly(integer 3) METhod(name min=1 max=1) lags(integer 999) TOLerance(real 0.00001) gmm /*
		  */ ATTrition ACF INIT(string) FSRESiduals(name min=1 max=1) EVALuator(string) TRANSlog OVERidentification]

	loc vv: di "version " string(min(max(10,c(stata_version)),14.0)) ", missing:" // we want a version 10 min, while a version 14.0 max (no version 14.2)
	loc depvar `varlist'
	marksample touse
	markout `touse' `free' `proxy' `state' `control' `endogenous' `id' `t'

	/// checks for same variables in state - free - proxy - controls
	loc chk1: list free & state
	loc chk2: list free & control
	loc chk3: list free & proxy
	loc chk4: list state & control
	loc chk5: list state & proxy
	loc chk6: list control & proxy
	forval i = 1/6{
		if "`chk`i''" != ""{
			di as error "Same variables in free, state, control or proxy."
			exit 198
		}
	}
	/// check whether there are more than one state AND more than one proxy
	loc pnum: word count `proxy'
	loc snum: word count `state'
	if `pnum' > 1 & `snum' > 1{
		di as error "Cannot specify multiple state AND multiple proxy variables"
		exit 198
	}
	/// check for unavailable choice of models
	if (!inlist("`method'", "op", "lp", "wrdg", "mr", "rob") & !mi("`method'")){
		di as error "Allowed prodest methods are op, lp, wrdg, mr or rob. Default is lp"
		exit 198
	}
	else if ("`method'" == "mr" & c(stata_version) < 14.2){
		di as error "MrEst only available with Stata version 14.2 or higher"
		exit 198
	}
	else if mi("`method'"){
		loc method "lp"
	}
	/// Check for unpractical value of polynomial approximation
	if (`poly' >= 7 | `poly' < 2){
		di as error "Polynomial degree must lie between 2 and 6"
		exit 198
	}
	/// Syntax check: is data xtset?
	if (mi("`id'") | mi("`t'")) {
	  capture xtset
	  if (_rc != 0) {
		 di as error "You must either xtset your data or specify both id() and t()"
		 error _rc
	  }
	  else {
		 loc id = r(panelvar)
		 loc t = r(timevar)
	  }
	}
	else {
	  qui xtset `id' `t'
	}
	/// Check for a valid confidence level
	if (`level' < 10 | `level' > 99) {
		di as error "confidence level must be between 10 and 99"
		error 198
	}
	/// Value added or gross output?
	if mi("`valueadded'"){
		loc model "grossoutput"
		loc proxyGO `proxy'
		loc lproxyGO l_gen`proxy'
	}
	else {
		loc model "valueadded"
	}
	/// Number of repetitions reasonable?
	if (`reps' < 2) & ("`method'" != "wrdg" & "`method'" != "mr"){
		di as error "reps() must be at least 2"
		exit 198
	}
	/// optimizer choice
	if (!mi("`optimizer'") & !inlist("`optimizer'", "nm", "nr", "dfp", "bfgs", "bhhh", "gn")){
		di as error "Allowed optimizers are nm, nr, dfp, bfgs and bhhh or gn for wrdg and mr. Default is nm or gn for wrdg and mr"
		exit 198
	}
	else if inlist("`method'", "wrdg", "mr") & inlist("`optimizer'", "nm", "bhhh"){
		di as error "`optimizer' is not allowed with `method' method. Optimizer switched to gn (default)"
		loc optimizer "gn"
	}
	else if mi("`optimizer'") & !inlist("`method'", "wrdg", "mr"){
		loc optimizer "nm"
	}
	else if mi("`optimizer'"){
		loc optimizer "gn"
	}
	/// number of lags in MR metodology
	if ("`method'" == "mr" & mi("`lags'")){
		di as error "Method 'mr' requires a specification for lags. Switched to 'all possible lags' (default)"
		loc lags = .
	}
	else if ("`method'" == "mr" & (`lags' < 1)) {
		di as error "Minimum lag is 1"
		exit 198
	}
	else if ("`method'" == "mr" & (`lags' == 999)) {
		loc lags = .
	}
	/// check ACF correction in woolridge and mr cases
	if ("`method'" == "wrdg" | "`method'" == "mr") & !mi("`acf'"){
		di as error "`method' does not support ACF correction. Estimation run on baseline model"
		loc acf ""
	}
	/// feasible values of tolerance
	if mi(`tolerance') | `tolerance' > 0.01{
		di as error "maximum value of tolerance is 0.01. Changed to default of value of e-5"
		loc tolerance = 0.00001
	}
	if ("`method'" == "op") {
		loc proxyGO ""
		loc lproxyGO ""
	}
	/// define conv_nrtol for wooldridge and mr - NOT in case of GN
	if ("`method'" == "wrdg" | "`method'" == "mr") & "`optimizer'" != "gn"{
		loc conv_nrtol "conv_nrtol(`tolerance')"
	}
	if ("`method'" == "wrdg" | "`method'" == "mr") & !mi("`init'"){
		loc init = subinstr("`init'",","," ",.)
		foreach var in `free' `state'{
			gettoken val init: init
			loc init_gmm `init_gmm' xb_`var' `val'
		}
		loc init_gmm from(`init_gmm')
	}
	if !mi("`fsresiduals'") & inlist("`method'", "wrdg", "mr", "rob"){
		di as error "fsresiduals is available with OP and LP methods only."
		exit 198
	}
	cap confirm var `fsresiduals'
	if !_rc{
		di as error "`fsresiduals' already exists"
		exit 198
	}
	/// define the evaluator type
	if mi("`evaluator'") & "`optimizer'" == "bhhh"{
		loc evaluator = "gf0"
	}
	else if mi("`evaluator'"){
		loc evaluator = "d0"
	}
	/// check the translog - only meaningful for ACF and Wooldridge methods
	if !mi("`translog'") & mi("`acf'") & "`method'" != "wrdg"{
		di as error "translog is available with Wooldridge or ACF-corrected models only"
		exit 198
	}
	else if !mi("`translog'") & !mi("`acf'"){
		loc transVars `free' `state' `proxyGO'
	}
	/// change the production function type
	if !mi("`translog'"){
		loc PFtype "translog"
	}
	else{
		loc PFtype "Cobb-Douglas"
	}
	if !mi("`overidentification'") & mi("`acf'") & "`method'" != "wrdg"{
		di as error "overidentification is meaningful for ACF or WRDG methods only. Ignoring the option."
	}
	if !mi("`gmm'") & "`method'" != "wrdg"{
		di as error "gmm works for wrdg only. Ignoring the option."
		loc gmm ""
	}
	if !mi("`acf'") & "`model'" == "grossoutput"{
		di as error "Using ACF correction with GO output does not ensure a correct parameter identification. See ACF (2015)."
	}

	if "`method'" == "op" loc strMethod "Olley-Pakes"
	if "`method'" == "lp" loc strMethod "Levinsohn-Petrin"
	if "`method'" == "wrdg" loc strMethod "Wooldridge"
	if "`method'" == "rob" loc strMethod "Wooldridge/Robinson"
	if "`method'" == "mr" loc strMethod "Mollisi-Rovigatti"

	loc colnum: word count `free' `state' `control' `proxyGO'

	/// initialize results matrices
	tempname firstb __b __V robV
	mat `__b' = J(`reps',`colnum',.)

	/// preserve the data before generating tons of variables
	preserve
		qui su `t' if `touse' == 1
		loc maxDate = `r(max)'

		/// directly keep only observations in IF and IN --> SAVE A TEMPORARY FILE
		tempfile temp
		qui keep if `touse' == 1
		keep `depvar' `free' `state' `proxy' `control' `id' `t' `touse'  `endogenous'
		/// generate an "exit" dummy variable equal to one for all firms not present in the last period of panel
		tempvar exit
		qui bys `id' (`t'): g `exit' = (_n == _N & `t' < `maxDate')
		qui save `temp', replace

		/// define all locals to run the command
		loc toLagVars `free' `state' `control' `proxyGO'

		foreach var of local toLagVars{
			qui g l_gen`var' = l.`var'
			loc laggedVars `laggedVars' l_gen`var'
		}
		foreach local in free state proxy control{
			loc `local'num: word count ``local''
			foreach var of local `local'{
				loc lag`local'Vars `lag`local'Vars' l_gen`var'
			}
		}

		loc instrumentVars `state' `lagfreeVars' `control' `lproxyGO'

		/// OP LP polyvars
		loc polyvars `state' `proxy'
		/// ACF requires free variables to be among the polynomial
		if !mi("`acf'"){
			loc polyvars `free' `polyvars'
		}
		loc varnum: word count `polyvars'
		loc controlnum: word count `control'
		loc tolagnum: word count `toLagVars'

		// poly-th degree polynomial
		loc n = 1
		foreach x of local polyvars{
			qui g var_`n' = `x'
			loc interactionvars `interactionvars' var_`n'
			loc ++n
		}

		forv i=1/`varnum'{
			forv j=`i'/`varnum'{
				qui g var_`i'_`j' = var_`i'*var_`j'
				loc interactionvars `interactionvars' var_`i'_`j'
				if `poly' > 2{
					forv z=`j'/`varnum'{
						qui g var_`i'_`j'_`z' = var_`i'*var_`j'*var_`z'
						loc interactionvars `interactionvars' var_`i'_`j'_`z'
						if `poly' > 3{
							forv g = `z'/`varnum'{
								qui g var_`i'_`j'_`z'_`g' = var_`i'*var_`j'*var_`z'*var_`g'
								loc interactionvars `interactionvars' var_`i'_`j'_`z'_`g'
								if `poly' > 4{
									forv v = `g'/`varnum'{
										qui g var_`i'_`j'_`z'_`g'_`v' = var_`i'*var_`j'*var_`z'*var_`g'*var_`v'
										loc interactionvars `interactionvars' var_`i'_`j'_`z'_`g'_`v'
										if `poly' > 5{
											forv s = `v'/`varnum'{
												qui g var_`i'_`j'_`z'_`g'_`v'_`s' = var_`i'*var_`j'*var_`z'*var_`g'*var_`v'*var_`s'
												loc interactionvars `interactionvars' var_`i'_`j'_`z'_`g'_`v'_`s'
											}
										}
									}
								}
							}
						}
					}
				}
			}
		}
		if inlist("`method'","wrdg", "mr", "rob"){
			/// generate GMM fit - for each state and free variables need to initiate the GMM
			foreach element in `free' `state' `proxyGO'{
				local gmmfit `gmmfit'-{`element'}*`element'
			}
			/// generate lag interactionvars
			foreach var of local interactionvars{
				cap g l_gen`var' = l.`var'
				loc lagInteractionvars `lagInteractionvars' l_gen`var'
			}
			if "`method'" == "wrdg"{ /* WRDG */
				if !mi("`overidentification'") | !mi("`gmm'") { // in case of overidentification, or in case the user specifies a preference for GMM, estimate a system GMM model, else go with the linear IV model
					if !mi("`overidentification'"){
						loc overidentification `lagfreeVars' `lagInteractionvars'
					}
					loc eq1counter = 0
					foreach var of varlist `interactionvars'{
						loc ++eq1counter
						loc eq1vars `eq1vars' -{xb`eq1counter'}*`var'
					}
					loc eq2counter = 0
					foreach var of varlist `lagInteractionvars'{
						loc ++eq2counter
						loc eq2vars `eq2vars' -{xb`eq2counter'}*`var'
					}
					`vv' qui gmm (eq1: `depvar' `gmmfit' `eq1vars' `wrdg_contr1' - {a0}) /* y - a0 - witB - xitJ - citL
						*/ (eq2: `depvar' `gmmfit' `eq2vars' `wrdg_contr2' - {a0} -{e0}), /* y - e0 - witB - xitJ - p(cit-1L) - ... - pg(cit-1L)^g
						*/ instruments(eq1: `free' `interactionvars' `control' `translogVars' `overidentification') /* Zit1 = (1,wit,xit,c0it,cit-1)
						*/ instruments(eq2: `state' `lagfreeVars' `lagInteractionvars' `control' `translogVars') /* Zit2 = (1,xit,wit-1,cit-1)
						*/ winitial(unadjusted, independent) technique(`optimizer') conv_maxiter(`maxiter') `conv_nrtol' `init_gmm'
					/// save locals for Hanses's J and p-value
					qui estat overid
					loc hans_j: di %3.2f `r(J)'
					loc hans_p: di %3.2f `r(J_p)'
				}
				else{ /* WRDG - plain */
					tempfile wrdg
					qui save `wrdg'
						qui reg `depvar' `state' `proxyGO' `interactionvars' `free' `lagfreeVars' `control' // just to take the number of observations used in the estimation
						loc realObs = `e(N)'
						qui expand 2, gen(cons2) // double the dataset in order to stack dependent and regressors
						foreach var of local interactionvars{
							qui replace `var' = l_gen`var' if cons2 == 1 // change interactionvars to lagged values
						}
						qui ivregress gmm `depvar' `state' `proxyGO' `interactionvars' `control' (`free' = `lagfreeVars') cons2, wmatrix(unadjusted) c // run the IV regression
						/// save locals for Hanses's J and its p-value
						loc hans_j: di %3.2f `e(J)'
						mat cV = colsof(e(V)) // find the number of instruments
						scalar cV = cV[1,1]
						loc jdf = cV - `e(rank)' + 1
						loc hans_p: di %3.2f chi2tail(`jdf', `hans_j') // this is the chi2 p-value of Hansen statistics
					qui use `wrdg', clear
					eret loc N `realObs' // post the real number of obseravtions used in estimation
				}
			}
			else if "`method'" == "mr"{ /* MrEst */
				if !mi("`overidentification'"){ // interactions are valid instruments, too
					loc overidentification  `lagInteractionvars'
				}
				/// god forgive me: in order to overcome the difference equation in system GMM we launch both equations in level and take the initial weighting matrix
				loc instnum = (`freenum' + `statenum')
				forv i = 1/`instnum'{
					loc intVars `intVars' var_`i'
				}
				loc interinstr: list local interactionvars - intVars
				qui gmm (`depvar' `gmmfit' - {xh: `interactionvars'} `wrdg_contr1' - {a0}), quickd instruments(1: `interinstr') /*
					*/ xtinstruments(`free' `state', l(0/`lags')) winitial(xt L) onestep conv_maxiter(1)
				qui mat W1 = e(W)
				qui gmm (`depvar' `gmmfit' - {xj: `lagInteractionvars'} `wrdg_contr2' - {a0} - {e0}), quickd /*
					*/ instruments(`state' `lagInteractionvars') xtinstruments(`free' `state', l(2/`lags'))/*
					*/ /*xtinstruments(`state', l(0/`lags'))*/ winitial(xt L) onestep conv_maxiter(1)
				qui mat W2 = e(W)
				mata W_hat =  st_matrix("W1"),J(rows(st_matrix("W1")),cols(st_matrix("W2")),0) \ J(rows(st_matrix("W2")),cols(st_matrix("W1")),0), st_matrix("W2")
				mata st_matrix("W_hat", W_hat)
				/* generate the model for both equations - same parameter for c_{i,t} and c_{i,t-1} in WRDG terms */
				loc eq1counter = 0
				foreach var of varlist `interactionvars' `control'{
					loc ++eq1counter
					loc eq1vars `eq1vars' "-{xb`eq1counter'}*`var'"
				}
				loc eq2counter = 0
				foreach var of varlist `lagInteractionvars' `control'{
					loc ++eq2counter
					loc eq2vars `eq2vars' "-{xb`eq2counter'}*`var'"
				}
				/// launch the gmm with the winitial built above
				`vv' qui gmm (1: `depvar' `gmmfit' `eq1vars' /*`wrdg_contr1'*/ - {a0}) /* y - a0 - witB - xitJ - citL
					*/ (2: `depvar' `gmmfit' `eq2vars' /*`wrdg_contr2'*/ - {a0} - {e0}),/* y - e0 - witB - xitJ - p(cit-1L) - ... - pg(cit-1L)^g, where p = g = 1
					*/ instruments(1: `interinstr' `overidentification' `control') xtinstruments(1: `free' `state', l(0/`lags')) /*
					*/ xtinstruments(2: `free' `state', l(2/`lags')) instruments(2: `state' `lagInteractionvars' `control') /*
					*/ onestep winitial("W_hat") nocommonesample quickd technique(`optimizer') conv_maxiter(`maxiter') `conv_nrtol' `init_gmm'
				/// save locals for Hanses's J and p-value
				qui estat overid
				loc hans_j: di %3.2f `r(J)'
				loc hans_p: di %3.2f `r(J_p)'
			}
			else{ /* Robinson / ACF */
				qui ivregress gmm `depvar' `state' `proxyGO' `control' `lagInteractionvars' (`free' = `lagfreeVars'), vce(cluster `id') /* this is working! */
				/// save locals for Hanses's J and its p-value
				loc hans_j: di %3.2f `e(J)'
				mat cV = colsof(e(V)) // find the number of instruments
				scalar cV = cV[1,1]
				loc jdf = cV - `e(rank)' + 1
				loc hans_p: di %3.2f chi2tail(`jdf', `hans_j') // this is the chi2 p-value of Hansen statistics
			}
			/// save elements for result posting
			loc nObs = `e(N)'
			loc numInstr1: word count `e(inst_1)'
			loc numInstr2: word count `e(inst_2)'
			mat `__b' = e(b)
			mat `__b' = `__b'[1...,1..`colnum']
			mat `__V' = e(V)
			mat `__V' = `__V'[1..`colnum',1..`colnum']
			/// save locals for Hanses's J and its p-value
			loc hans_j: di %3.2f `e(J)'
			mat cV = colsof(e(V)) // find the number of instruments
			scalar cV = cV[1,1]
			loc jdf = cV - `e(rank)' + 1
			loc hans_p: di %3.2f chi2tail(`jdf', `hans_j') // this is the chi2 p-value of Hansen statistics
			continue, break
		}
		else{ /* if it's not WRDG or MrEst */
			if !mi("`attrition'"){
				foreach var in `interactionvars'{
					qui g l_gen`var' = l.`var'
					loc lagInterVars `lagInterVars' l_gen`var'
				}
				qui cap logit `exit' `lagInterVars'
				if _rc == 0{
					qui predict Pr_hat if e(sample), pr
				}
				else{
					di as error "No ID exits the sample. Running the estimation with no attrition"
					loc attrition ""
					qui g Pr_hat = .
				}
			}
			else{
				qui g Pr_hat = .
			}
			/// in case of ACF we don't want free variables to appear twice in the regression
			if !mi("`acf'"){
				forv i = 1/`freenum'{
					loc freeVars `freeVars' var_`i'
				}
				loc interactionvars: list local interactionvars - freeVars
				/// generate the needed variables for translog: a local with all power of 2 + interactions and instruments - lag of free/proxy interacted with state
				if !mi("`translog'"){
					loc transNum: word count `transVars'
					forv i=1/`transNum'{
						forv j=`i'/`transNum'{
							loc ivar `: word `i' of `transVars''
							loc jvar `: word `j' of `transVars''
							loc tvarname "`ivar'X`jvar'"
							loc transNames `transNames' `tvarname'
							loc interactionTransVars `interactionTransVars' var_`i'_`j'
							qui g lagTrans_`i'_`j' = l.var_`i'_`j'
							loc lagInteractionTransVars `lagInteractionTransVars' lagTrans_`i'_`j'
						}
					}
					foreach fvar in `free' `proxyGO'{
						tempvar d`fvar'
						qui g `d`fvar'' = `fvar'^2
						qui g lagInstr`fvar' = l.`d`fvar''
						loc instrumentTransVars `instrumentTransVars' lagInstr`fvar'
						foreach svar in `state'{
							qui g instr_`fvar'`svar' = lagInstr`fvar'*`svar'
							loc instrumentTransVars `instrumentTransVars' instr_`fvar'`svar'
						}
					}
					foreach svar in `state'{
						tempvar d`svar'
						qui g `d`svar'' = `svar'^2
						qui g lagInstr`svar' = l.`d`svar''
						loc instrumentTransVars `instrumentTransVars' lagInstr`svar'
					}
					loc colnum: word count `free' `state' `control' `proxyGO' `interactionTransVars'
					mat `__b' = J(`reps',`colnum',.)
				}
			}
			loc regvars `free' `control' `interactionvars'
			loc firstRegNum: word count `regvars'
			loc regNum: word count `free' `state' `control' `proxyGO' `transVars'

			/// first stage
			qui _xt, trequired
			qui reg `depvar' `regvars'
			/// generating "freeFit" as the fit of free variables to be subtracted to Y
			tempvar freeFit
			qui g `freeFit' = 0
			foreach var in `free'{
				scalar b_`var' = _b[`var']
				qui replace `freeFit' = `freeFit' + (b_`var'*`var')
			}
			mat `firstb' = e(b)
			mat `robV' = e(V)
			qui predict phihat if e(sample), xb
			qui g phihat_lag = l.phihat

			/// retrieve starting points through an OLS estimation --> first round only, not during bootstrap
			qui reg `depvar' `toLagVars' `interactionTransVars' if `touse' == 1
			if !mi("`init'"){
				mat ols_s = `init' // THIS PART IS JUST A TRYOUT IN ORDER TO MAKE THE COMMAND WORK FOR OUR PURPOSES
			}
			else{
				mat tmp = e(b)
				mata: st_matrix("ols_s",st_matrix("tmp"):+ rnormal(1,1,0,.01)) // we add some "noise" to OLS results in order not to block the optimizer
			}
			/// save the first stage results
			if !mi("`fsresiduals'"){
				tempfile fsres
				tempvar FSfit
				mat score `FSfit' = `firstb'
				qui g `fsresiduals' = `depvar' - `FSfit'
				qui save `fsres'
			}
			qui g res = 0
			if mi("`acf'"){ /* OP and LP (non-corrected) second stage */
				/// here we generate a tempvar with the fitted value of all the free variables
				qui replace phihat = phihat - `freeFit'
				qui replace res = `depvar' - `freeFit'
				qui replace phihat_lag = l.phihat
				mat init = ols_s[1...,(`freenum'+1)..`regNum']'
				loc toLagVars `state' `control' `proxyGO'
				loc laggedVars: list local laggedVars - lagfreeVars
				loc instrumentVars `state' `control' `lproxyGO'
				/// here we launch the mata routine for OP or LP
				foreach var of varlist `toLagVars' `laggedVars' /*`lagfreeVars'*/ phihat_lag phihat{
					qui drop if mi(`var')
				}
				/// the routine cannot fail the first estimation - otherwise it would fail with the point estimates. We capture the boot repetitions
				qui mata: opt_mata(st_matrix("init"),&foplp(),"`optimizer'","phihat","phihat_lag",/*
					*/"`toLagVars'","`laggedVars'","`touse'",`maxiter',`tolerance',"`evaluator'","Pr_hat","res","`instrumentVars'","`endogenous'")
				mat `__b'[1,1] = `firstb'[1...,1..(`freenum')],r(betas)
			}
			else{ /* ACF second stage */
				if !mi("`overidentification'"){ // lag intereactions are valid instruments for the first equation, too
					loc overidentification `lagInteractionvars'
				}
				loc toLagVars `free' `state' `control' `proxyGO' `interactionTransVars'
				loc laggedVars `lagfreeVars' `lagstateVars' `lagcontrolVars' `lproxyGO' `lagInteractionTransVars'
				loc instrumentVars `lagfreeVars' `state' `control' `lproxyGO' `instrumentTransVars' `overidentification'
				/// here we launch the mata routine for ACF
				foreach var of varlist `laggedVars' `lagfreeVars' phihat_lag{
					qui drop if mi(`var')
				}
				loc betaNum: word count `free' `state' `proxyGO' `control' `interactionTransVars'
				mat init = ols_s[1...,1..(`betaNum')]'
				qui mata: opt_mata(st_matrix("init"),&facf(),"`optimizer'","phihat","phihat_lag",/*
					*/ "`toLagVars'","`laggedVars'","`touse'",`maxiter',`tolerance',"`evaluator'","Pr_hat","res","`instrumentVars'","`endogenous'")
				mat `__b'[1,1] = r(betas)
			}
		}

		if !inlist("`method'","wrdg","mr","rob"){
			clear
			/// generate the varCovar matrix for the bootstrapped estimates
			qui svmat `__b'
			mat `__b' = `__b'[1,1...]
		}

	restore

	/// merge the first stage residuals
	if !mi("`fsresiduals'"){
		qui merge 1:1 `id' `t' using `fsres', nogen keepusing(`fsresiduals')
	}

	mat coleq `__b' = ""
	mat colnames `__b' = `free' `state' `control' `proxyGO' `transNames'

	/// Display results - ereturn
	eret clear

	eret post `__b'

	if !mi("`acf'"){
		loc correction "ACF corrected"
	}
	if !mi("`fsresiduals'"){
		eret loc FSres "`fsresiduals'"
	}

	eret loc cmd "prodest"
	eret loc depvar "`depvar'"
	eret loc free "`free'"
	eret loc state "`state'"
	eret loc proxy "`proxy'"
	eret loc controls "`control'"
	eret loc endogenous "`endogenous'"
	eret loc method "`strMethod'"
	eret loc model "`model'"
	eret loc technique "`optimizer'"
	eret loc idvar "`id'"
	eret loc timevar "`t'"
	eret loc correction "`correction'"
	eret loc predict "prodest_p_m"
	eret loc PFtype "`PFtype'"
	eret loc gmm "`gmm'"

end program


//*---------------------------------------------------------------------*/

/// defining mata routines for optimization
capture mata mata drop opt_mata()
capture mata mata drop facf()
capture mata mata drop foplp()


*version 7
mata:

/*---------------------------------------------------------------------*/

	void foplp(todo,betas,X,lX,PHI,LPHI,RES,Z,PR_HAT,ENDO,crit,g,H)
	{
		OMEGA = PHI-X*betas'
		OMEGA_lag = LPHI-lX*betas'
		OMEGA_lag2 = OMEGA_lag:*OMEGA_lag
		OMEGA_lag3 = OMEGA_lag2:*OMEGA_lag
		/* IF clause in order to see whether we have to use the "exit" variable */
		if (!missing(PR_HAT)){
			PR_HAT2 = PR_HAT:*PR_HAT
			PR_HAT3 = PR_HAT2:*PR_HAT
			OMEGA_lag_pol = (J(rows(PHI),1,1),OMEGA_lag,OMEGA_lag2,OMEGA_lag3,PR_HAT,PR_HAT2,PR_HAT3,PR_HAT:*OMEGA_lag,PR_HAT2:*OMEGA_lag,PR_HAT:*OMEGA_lag2,ENDO)
		}
		else{
			OMEGA_lag_pol = (J(rows(PHI),1,1),OMEGA_lag,OMEGA_lag2,OMEGA_lag3,ENDO)
		}
		g_b = invsym(OMEGA_lag_pol'OMEGA_lag_pol)*OMEGA_lag_pol'OMEGA
		XI = RES-X*betas'-OMEGA_lag_pol*g_b
		crit = (XI)'*(XI)
	}

/*---------------------------------------------------------------------*/

	void facf(todo,betas,X,lX,PHI,LPHI,RES,Z,PR_HAT,ENDO,crit,g,H)
	{
		W = invsym(Z'Z)/(rows(Z))
		OMEGA = PHI-X*betas'
		OMEGA_lag = LPHI-lX*betas'
		OMEGA_lag2 = OMEGA_lag:*OMEGA_lag
		OMEGA_lag3 = OMEGA_lag2:*OMEGA_lag
		/* IF clause in order to see whether we have to use the "exit" variable */
		if (!missing(PR_HAT)){
			PR_HAT2 = PR_HAT:*PR_HAT
			PR_HAT3 = PR_HAT2:*PR_HAT
			OMEGA_lag_pol = (J(rows(PHI),1,1),OMEGA_lag,OMEGA_lag2,OMEGA_lag3,PR_HAT,PR_HAT2,PR_HAT3,PR_HAT:*OMEGA_lag,PR_HAT2:*OMEGA_lag,PR_HAT:*OMEGA_lag2,ENDO)
		}
		else{
			OMEGA_lag_pol = (J(rows(PHI),1,1),OMEGA_lag,OMEGA_lag2,OMEGA_lag3,ENDO)
		}
		g_b = invsym(OMEGA_lag_pol'OMEGA_lag_pol)*OMEGA_lag_pol'OMEGA
		XI = OMEGA-OMEGA_lag_pol*g_b
		crit = (Z'XI)'*W*(Z'XI)
	}

/*---------------------------------------------------------------------*/

void opt_mata(init, f, opt, phi, lphi, tolag, lagged, touse, maxiter, tol, eval, | Pr_hat, res, instr, endogenous)
	{
		st_view(RES=.,.,st_tsrevar(tokens(res)), touse)
		st_view(PHI=.,.,st_tsrevar(tokens(phi)), touse)
		st_view(LPHI=.,.,st_tsrevar(tokens(lphi)), touse)
		st_view(Z=.,.,st_tsrevar(tokens(instr)), touse)
		st_view(X=.,.,st_tsrevar(tokens(tolag)), touse)
		st_view(lX=.,.,st_tsrevar(tokens(lagged)), touse)
		st_view(PR_HAT=.,.,st_tsrevar(tokens(Pr_hat)), touse)
		st_view(ENDO=.,.,st_tsrevar(tokens(endogenous)), touse)
		S = optimize_init()
		optimize_init_argument(S, 1, X)
		optimize_init_argument(S, 2, lX)
		optimize_init_argument(S, 3, PHI)
		optimize_init_argument(S, 4, LPHI)
		optimize_init_argument(S, 5, RES)
		optimize_init_argument(S, 6, Z)
		optimize_init_argument(S, 7, PR_HAT)
		optimize_init_argument(S, 8, ENDO)
		optimize_init_evaluator(S, f)
		/*optimize_init_evaluatortype(S, "d0")*/
		optimize_init_evaluatortype(S, eval)
		optimize_init_conv_maxiter(S, maxiter)
		optimize_init_conv_nrtol(S, tol)
		/*optimize_init_evaluatortype(S, "gf0")*/
		optimize_init_technique(S, opt)
		optimize_init_nmsimplexdeltas(S, 0.00001)
		optimize_init_which(S,"min")
		optimize_init_params(S,init')
		p = optimize(S)
		st_matrix("r(betas)", p)

		/* tryout to find Newey-type Cov estimator
		st_matrix("r(gradient)", optimize_result_gradient(S))
		st_matrix("r(score)", optimize_result_scores(S))
		Hl = optimize_result_gradient(S)
		V = optimize_result_V(S)
		Vhh = invsym(Hl)*V*invsym(Hl)'
		st_matrix("r(Vhh)", Vhh)*/
	}

/*---------------------------------------------------------------------*/

end