quietly{
set mem 10m

* Load Data
use GMdata

log using PS1, replace text

* QUESTION 1
noisily display("Question 1")
* Get number of periods firm is in the sample
bysort index: egen perinsamp = count(index)
iis index
tis yr
* Get means and standard deviations
noisily sum l*
noisily sum l* if perinsamp==4
* If in the sample at least 2 periods
noisily sum l* if perinsamp>=2

* QUESTION 2
noisily display("Question 2")
*Generate computer industry indicator
gen d357 = (sic==357)
* i) regressions:
* Total
xi:regress ldsal lemp ldnpt ldrst i.yr*i.d357 if perinsamp==4, robust
estimates store ols1

* Between
xi:xtreg ldsal lemp ldnpt ldrst i.yr*i.d357 if perinsamp==4, be 
estimates store olsbe

* Within
xi:xtreg ldsal lemp ldnpt ldrst i.yr*i.d357 if perinsamp==4, fe robust
estimates store panelwithin

* Random
xi:xtreg ldsal lemp ldnpt ldrst i.yr*i.d357 if perinsamp==4, robust
estimates store panelrandom

noisily estimates table ols1 olsbe panelwithin panelrandom, se

* ii) Hausman Test:
noisily hausman panelwithin panelrandom


* QUESTION 3
noisily display("Question 3")

* Total 
xi:regress ldsal lemp ldnpt ldrst i.yr*i.d357  
estimates store olsub1

* First Difference
xi:regress S5.ldsal S5.lemp S5.ldnpt S5.ldrst i.yr*i.d357 
estimates store olsub2

* At least 2 periods Total
xi:regress ldsal lemp ldnpt ldrst i.yr*i.d357 if perinsamp>=2
estimates store olsub3

noisily estimates table olsub1 olsub2 olsub3, se

* ii) Survival Regression
sort index yr
by index: generate death=index[_n+1]==.
replace death = . if yr==88 
generate life = 1-death
probit life ldnpt ldrst ldinv
estimates store theprobit

noisily display("Exit Regression")
noisily estimates replay theprobit
predict lifeprob

* Inverse Mills
generate imrlife=-normalden(lifeprob)/normal(lifeprob)

* At least 2 periods Total
xi:regress ldsal lemp ldnpt ldrst i.yr*i.d357 imrlife if perinsamp>=2
estimates store heck1

tsset index yr
* First Difference at least 2 periods
xi:regress S5.ldsal S5.lemp S5.ldnpt S5.ldrst i.yr*i.d357 imrlife  if perinsamp>=2
estimates store heck2
noisily display("Heckman Type estimates")
noisily estimates table heck1 heck2, se

/* Now there is also a simpler way of doing this in Stata given the HECKMAN command, but it is not
really in the spirit of the question */
xi:heckman ldsal lemp ldnpt ldrst i.yr*i.d357 if perinsamp>=2, select(death=ldnpt ldrst ldinv)
estimates store heck1
xi:heckman S5.ldsal S5.lemp S5.ldnpt S5.ldrst i.yr*i.d357 if perinsamp>=2, select(death=ldnpt ldrst ldinv)
estimates store heck2
noisily display("One line estimates of Heckman Model")
noisily estimates table heck1 heck2, se


* QUESTION 4
noisily display("Question 4")
tsset index yr

* part (i)
local bests = 100
local bestrho = 0
forvalues rho= 0.2 (0.001) .8 {
gen qdy = ldsal - `rho'*S5.ldsal
gen qdl = lemp - `rho'*S5.lemp
gen qdk = ldnpt - `rho'*S5.ldnpt
gen qdr = ldrst - `rho'*S5.ldrst

xi: ivreg qdy (qdl qdk qdr = S5.lemp S5.ldnpt S5.ldrst)i.yr*i.d357
predict resid, residuals
gen x1 = resid*S5.ldsal 
summ x1
local s1 = abs(r(N)*r(mean))

if `s1' < `bests' {
local bests = `s1'
local bestrho = `rho'
}
drop qdy qdl qdk qdr resid x1
}
noisily macro list _bestrho _bests

gen qdy = ldsal - `bestrho'*S5.ldsal
gen qdl = lemp - `bestrho'*S5.lemp
gen qdk = ldnpt - `bestrho'*S5.ldnpt
gen qdr = ldrst - `bestrho'*S5.ldrst

noisily xi: ivreg qdy (qdl qdk qdr = S5.lemp S5.ldnpt S5.ldrst)i.yr*i.d357
drop qdy qdl qdk qdr


* (ii)
local bests = 100
local bestrho = 0
forvalues rho= 0.2 (0.001) .8 {
gen qdy = (ldsal - `rho'*S5.ldsal) - (S5.ldsal - `rho'*S10.ldsal) 
gen qdl = (lemp - `rho'*S5.lemp) - (S5.lemp - `rho'*S10.lemp)
gen qdk = (ldnpt - `rho'*S5.ldnpt) - (S5.ldnpt - `rho'*S10.ldnpt)
gen qdr = (ldrst - `rho'*S5.ldrst) - (S5.ldrst - `rho'*S10.ldrst)

xi: ivreg qdy (qdl qdk qdr = S5.lemp S5.ldnpt S5.ldrst)i.yr*i.d357
predict resid, residuals
gen x1 = resid*S10.ldsal 
summ x1
local s1 = abs(r(N)*r(mean))

if `s1' < `bests' {
local bests = `s1'
local bestrho = `rho'
}
drop qdy qdl qdk qdr resid x1
}
noisily macro list _bestrho _bests

gen qdy = (ldsal - `bestrho'*S5.ldsal) - (S5.ldsal - `bestrho'*S10.ldsal) 
gen qdl = (lemp - `bestrho'*S5.lemp) - (S5.lemp - `bestrho'*S10.lemp)
gen qdk = (ldnpt - `bestrho'*S5.ldnpt) - (S5.ldnpt - `bestrho'*S10.ldnpt)
gen qdr = (ldrst - `bestrho'*S5.ldrst) - (S5.ldrst - `bestrho'*S10.ldrst)

noisily xi: ivreg qdy (qdl qdk qdr = S5.lemp S5.ldnpt S5.ldrst)i.yr*i.d357
drop qdy qdl qdk qdr



* QUESTION 5: OLLEY-PAKES
noisily display("Question 5")
* i)

gen ldnptsq = ldnpt*ldnpt
gen ldrstsq = ldrst*ldrst
gen ldinvsq = ldinv*ldinv
gen ldnptldrst=ldnpt*ldrst
gen ldnptldinv=ldnpt*ldinv
gen ldrstldinv=ldrst*ldinv
noisily display("Regression on year dummies and polynomial")
noisily xi:reg ldsal lemp i.yr*i.d357 ldnpt* ldrst* ldinv* 


* Generate new dependent variable 
tsset index yr
gen newy = F5.ldsal-_b[lemp]*F5.lemp-_b[_Iyr_78]*F5._Iyr_78-_b[_Iyr_83]*F5._Iyr_83-_b[_Iyr_88]*F5._Iyr_88-_b[_Id357_1]*F5._Id357_1-_b[_IyrXd35_78_1]*F5._IyrXd35_78_1-_b[_IyrXd35_83_1]*F5._IyrXd35_83_1-_b[_IyrXd35_88_1]*F5._IyrXd35_88_1 
gen polynew = _b[ldnpt]*ldnpt+_b[ldnptsq]*ldnptsq+_b[ldrst]*ldrst+_b[ldrstsq]*ldrstsq+_b[ldinv]*ldinv+_b[ldinvsq]*ldinvsq+_b[ldnptldrst]*ldnptldrst+_b[ldnptldinv]*ldnptldinv+_b[ldrstldinv]*ldrstldinv 

* ii) Non-linear Least Squares

nl (newy = {b0}+{b2}*F5.ldnpt+{b3}*F5.ldrst+{bh}*(polynew-{b2}*ldnpt-{b3}*ldrst)+{bhsq}*(polynew-{b2}*ldnpt-{b3}*ldrst)^2), 
estimates store olleypakes
noisily estimates replay olleypakes

*iii) Propensity Scores
gen lifeprobsq = lifeprob*lifeprob

noisily display("Regression on Phat and Phatsquared")
noisily reg newy ldnpt ldrst lifeprob lifeprobsq


*iv) Mixed polynomial propensity score
nl (newy = {b0}+{b2}*F5.ldnpt+{b3}*F5.ldrst+{bh}*(polynew-{b2}*ldnpt-{b3}*ldrst)+{bhsq}*(polynew-{b2}*ldnpt-{b3}*ldrst)^2+{bP}*lifeprob+{bPsq}*lifeprobsq+{bhP}*lifeprob*(polynew-{b2}*ldnpt-{b3}*ldrst))
estimates store opp
noisily estimates replay opp

} 
log close
display("end")