Lorenz Attractor Simulation

# This is a brief exploration of Chaos theory.  The Lorenz attractor is an equation used to model convection.  It is defined by the system of equations:

# dx/dt = sigam(y-x)

# dy/dt = x(rho-z)-y

# dz/dt = xy-beta*z

# With sigma, rho, and beta representing model parameters.

# We can easily approximate this system by a series of discreet time steps

# First set the initial values

y = 5

x = 5

z = 5

# Now let's set the parameters

sigma = 10

rho = 28

beta = 8/3

for (i in 1:9999) {

  x[i+1] = x[i] + sigma*(y[i]-x[i])/200

  y[i+1] = y[i] + (x[i]*(rho-z[i])-y[i])/200

  z[i+1] = z[i] + (x[i]*y[i]-beta*z[i])/200

}

plot(x[!is.na(x)],y[!is.na(x)], type="n")

for(i in 1:(length(x)-1)) 

{ 

lines (x[i:(i+1)], y[i:(i+1)], col = rainbow(length(x))[i]) 

} 

Law of Iterative Expectations/Law of Total Expectations

* The law of iterative expectations (LIE) is an extremely useful rule that can be frequently and usefully employed in many econometric proofs.

* This post will explore some simulated results in an attempt to strengthen our intution.

* LIE is the statement E(Y)=E(E(Y|X))

* For a proof using sums (http://en.wikipedia.org/wiki/Law_of_total_expectation)

* A proof using pdf probability rules can be found at (http://econ.wikidot.com/conditionaldistributions)

* Example 1: OLS

* In OLS (given the linearity assumption and the zero conditional mean assumption) E(Y|X) = XB

* Thus using LIE: E(Y) = E(E(Y|X)) -> in the sample analogue mean(Y) = mean(XBhat) = mean(Yhat)

* Let's see this in action

clear
set obs 1000

gen x = rnormal() + 5
gen u = rnormal()

gen y = x + u*10

reg y x

predict yhat
sum y yhat

* Example 2: Imagine you are trying to calculate the reccomended cooking times for a frozen food that you have developed that is best for the most people.

* You know that for each 1000 feet above sea level you need to reduce the temperature by 5% (this is entirely fictional).  You have experimented and found the ideal cooking temperature at sea level is 350. The result is E(Y|x)=350*(.95)^x where x is 1000 feet.

* Now you want to calculate a single cooking tempurature reccomendation that gives the ideal tempurature for an entire population.

* You know 30% of consumers live at sea level, 20% at 500 feet, 20% at 1000, 10% at 2000, 10% at 3500, and 10% at 5000.

* Thus expected best temperature reccomendation for a randomly drawn person is:

di .3*350*(.95)^0 + .2*350*(.95)^.5 + .2*350*(.95)^1 + .1*350*(.95)^2 + .1*350*(.95)^3.5 + .1*350*(.95)^5

* Thus E(Y)=E(E(Y|X))

* Note that this is a different temperature than taking the temperature at the average elevation.

di "Average elevation = " .3*0 + .2*500 + .2*1000 + .1*2000 + .1*3500 + .1*5000

di "Temperature reccomendation at average elevation = " 350*(.95)^((.3*0+.2*500+.2*1000+.1*2000+.1*3500+.1*5000)/1000)

Inverse Probability Wieghting to Correct for Sample Selection/Missing Data

* Imagine that you have some data set that is missing some of the variables of interest but you have a complete set of explanatory variables.  You might be concerned that the selection from the sample is correlated or is causing correlation in errors with your explanatory variable of interests thus creating potential bias.

* Imagine that you are interested in estimating if obedience school for dogs has the potential to reduce their risk of biting people.  As the y variable you have self-reported (by owners) number of bites.  As the explanatory variable you have breed aggressiveness and an indicator if the dog went to obedience school.

* Imagine also that you have information on the aggressiveness of the owners of the dogs which is correlated with the error in estimating the number of bites.  It is also an explanatory variable of selection.

clear
set obs 100000

gen n_classes = rpoisson(1)
gen breed_agg = rnormal()
gen owner_agg = rnormal()

* First lets calculate selection
gen p = normal(.5 -.5*n_classes + .5*owner_agg)
gen s = rbinomial(1,p)
gen u = rnormal()+owner_agg

* First let's assume there is no selection
gen bites = 2 + breed_agg - n_classes + 3*u

reg bites breed_agg n_classes
* We can see in this case the estimates look good (absent of selection)

replace bites = . if s == 0

reg bites breed_agg n_classes
* Now, the effects of classes seem greatly diminished in our observables because of the correlation between selection and the error component resulting from ownernship aggressiveness.

* There are two ways I can think of generating an unbiased estimates.

* We must do this by removing the correlation with selection and the error.

reg bites breed_agg n_classes owner_agg
* Is the easiest way to do this.  However, this post is about inverse probability weighting.  So that is what we will do.

probit s n_classes owner_agg

predict shat
* We want to estimate probability of selection from observables

gen ishat = 1/shat
* Then first the inverse of it to use in the pweight command

reg bites breed_agg n_classes [pweight=ishat]
* We can see this estimate is working well though the previous regression was generally better.

* This post deals with inverse probability weighting in simple OLS.  A future post will address inverse probability weighting in M-estimation: http://ideas.repec.org/p/ifs/cemmap/11-02.html

The DELTA method

# The delta method is an extremely useful tool for estimating the standard errors of non-linear models.

# This post will use as a reference David Patterson of the University of Montana's post on the delta method http://www.math.umt.edu/patterson/549/Delta.pdf

# The delta method states that the standard error of a distribution f(x) is approximated by f'(Mx) var(x) f(Mx) where f(x) is the gradiant of F(x) with respect to x (first derivative), and Mx=mean(x).

# Let's first generate a base variable x~N with mean 5 and standard deviation 3.
x=rnorm(10000)*3+5

# Example 1: F(x) = x^2 -> f(x) = 2x ->
# thus var(F(x))~2*mean(x)*var(x)*2*mean(x) ->
# var(F(x)) ~
  2*5*9*2*5

# Let's see how well our approximation compares with actual data
x2 = x^2
var(x2)

# We can see that our estimate (900) is pretty close to the sample variance (1050ish)

# Using our sampling distribution
2*mean(x)*var(x)*2*mean(x)

# Example 2: F(x) = 1/x -> f(x) = (-1)*x^(-2)
# var(F(x)) ~ (-1)*mean(x)^(-2) * var(x) * (-1)*mean(x)^(-2)
# = mean(x)^(-4) * var(x)
mean(x)^(-4) * var(x)

# Let's see how well our approximation compares with actual data
xi = 1/x

var(xi)
# In this case, not very well.  This is probably because the distribution of x approaches the nondiferentiable point 0 for some xs making a first order approximation fail badly .

# Let us try with another variable
z = rnorm(10000)+20

mean(z)^(-4) * var(z)

zi = 1/z

var(zi)
# We can see now the delta method is working well.

# Example 3: F(z) = ln(z) -> f(z) = 1/z
# var(F(z)) ~ (1/mean(z))*var(z)*(1/mean(z))
# = (1/mean(z))^2*var(z)
(1/mean(z))^2*var(z)

lnz = log(z)

var(lnz)
# Once again, looking good.

# Now let's imagine an estimator.  First off we need to recognize that estimators are random variables.

# Bhat = (X'X)^(-1)X'Y = (X'X)^(-1)X'(XB + e) = (X'X)^(-1)X'XB + (X'X)^(-1)X'e
#      = B + (X'X)^(-1)X'e
# E(Bhat|x) = E(B|x) + E[(X'X)^-1 X'e|x] -> assume (e|x)=0 ->
# E(Bhat|x) = B + (X'X)^-1 X'E[e|x] = B
# var(Bhat|x) = var(B + (X'X)^(-1)X'e|x) = var((X'X)^(-1)X'e|x)
# = (X'X)^(-1)X'var(e|x)((X'X)^(-1)X')' = (X'X)^(-1)X'var(e|x)X(X'X)^(-1)
# = (X'X)^(-1)X' *I X(X'X)^(-1)
# = sigma2*(X'X)^(-1)X'X(X'X)^(-1) = sigma2*(X'X)^(-1)

e = rnorm(1000)
y = 5*x + e*100

# Let's estimate OLS without residuals
ols=lm(y~x-1)

summary(ols)

uhat = residuals(ols)

# var(Bhat|x) ~ sigma2(X'X)^(-1) estimate
var(uhat)*sum(x*x)^(-1)

# standard error
(var(uhat)*sum(x*x)^(-1))^.5

Linear models with heteroskedastic errors

* Imagine we have a sample of hotels and prices of rooms in those hotels
clear
set obs 200
set seed 101
* This sets the random seed at 101.  Thus every time this simulation is run it will product identical results.

gen star = ceil(runiform()*7)/2+.5
  label var star "Number of stars of hotel"

gen id = _n

gen het = abs(rnormal())
  label var het "Unobserved heterogeneity."

* Generate a the number of reservations observed per hotel.
gen num_reservations = rpoisson(150)
expand num_reservations

gen seasonality = abs(rnormal())
  label var seasonality "Seasonal demand"

gen v = abs(rnormal())
  label var v "Unobserved variance in the variance term"

gen u = rnormal()*(5+het*15+7.5*star+22.5*seasonality+v*10)
  label var u "Error term"

gen p = 175 + 20*star + 15*seasonality + u + het

sum p
* There are some p values which are less than 0 but we can think of those as special deals, coupons, refunds, or other situations that might result in the effective price being less than 0 dollars.
* If we were to eliminate the less than 0 prices then we would in be enforcing left censoring which is a different problem.  See "tobit".  This blog has several posts touching on the use of the tobit.

* Estimate the price through direct OLS
reg p star seasonality

* Set the panel level observation
xtset id

* Though we have panel data we cannot effectively use fixed effect or random effects approaches to identify the price effect having one more star has on prices.
xtreg p star seasonality, fe
xtreg p star seasonality, re

* Now let's attempt a two step method to more efficient identify the error and reestimate the OLS.

reg p star seasonality
* The OLS regression looks pretty good.  Let's see if we can improve on it.  Note, the 95% confidence interval did not capture the true coefficient of 20 but that is not necessarily a problem.

reg p star seasonality, robust
reg p star seasonality, cluster(id)
* Using robust and cluster robust estimates of the standard error does not change the variance so much that the 95% confidence interval encloses the 20.

predict uhat, resid

gen uhat_abs = abs(uhat)
  label var uhat_abs "Abs of OLS residual"

two  (scatter uhat_abs seasonality) (scatter uhat_abs star), ///
              legend(label(1 "Seasonality") label(2 "Stars"))

* Since E(u)==0 var(u) is equal to u squared Var(u)=(u-E(u))^2
* Likewise, u^2 is approximated by u^2
* Similarly sd(u) should be approximated by (u^2)^.5=abs(u)
reg uhat_abs star seasonality

predict uhat_abs_hat, xb
  * I might be doing this wrong.  I am trying aweights which say something about weighting by the inverse of the variance of the observation.)

gen uhat2 = 1/uhat_abs_hat^2

reg p star seasonality
* Let's see how the unweighted estimate performs

reg p star seasonality [aweight = uhat2]
* Using variance weights does not seem to improve the estimate

* Alternatively we can use the MLE estimator allowing the conditional standard deviation of the error as well as the conditional mean to vary linearly.

cap program drop mle_ols
program mle_ols
  args log_like xb sigma
  qui replace `log_like' = ln(normalden($ML_y1-`xb',0,`sigma'))
end

ml model lf mle_ols (price: p = star seasonality) (sigma: star seasonality)
ml maximize

* Using the MLE estimator we seem to have gained precision in the 3rd decimal place of the coefficients.  The 95% CI still does not enclose the 20 but it is closer.  Still, this is not indicative of a problem.  Perhaps if we simulated this 1000 times and substantially more than 50 times the CI did not enclose the 20 then we might be worried.

A brief look at fe vs re

global num_reps = 100
* Set the number of repetitions

quietly forv i = 1(1)$num_reps {
* Loop through the simulation 100 times

gl noi
if `i'==1 gl noi noi
  * Only display regression results on the firt loop

clear
set obs 1000

gen id = _n

gen het = rnormal()
  label var het "Individual heterogeneity"
gen xauto = rnormal()
  label var xauto "Autocorrelation in x"
gen uauto = rnormal()
  label var uauto "Autocorrelatoin in u"

expand 5
  * Create 5 observations per individual

gen x = rnormal() + .5*xauto
  * Generate an x that is part auto correlated and part unique draws
 
gen y = 5*x + rnormal()*50 + uauto*25
  * Generate a y in which error is also autocorrelated by individual

bysort id: gen t=_n
  * Create a time variable

xtset id t
  * Assign panel data indicators

$noi xtreg y x, fe
  * Run a fixed effect regression
gl fe`i' =  _b[x]
  * Save the results in the global fe#
 

$noi xtreg y x, re
  * Run a random effect regression
gl re`i' =  _b[x]
  * Save the results in the global re#

noi di `i'
  * Display repetition number
}

* Clear the old memory
clear
set obs $num_reps

* Create empty variable holders.
gen re=.
gen fe=.

* Loop through the repetititions saving each to a variable
forv i=1(1)$num_reps {
  replace re = ${re`i'} if _n==`i'
  replace fe = ${fe`i'} if _n==`i'
}

sum re fe

Proof IV is same as 2SLS

This is a standard proof.  Normally I am the last one to attempt proofs but I wanted to try out $ \LaTeX{} $ in blogger.






































LaTeX Script:
\documentclass[10pt]{report}
\usepackage[utf8]{inputenc}
\usepackage{amsmath}
\begin{document}

  % This is a comment; it is not shown in the final output.
  % The following shows a little of the typesetting power of LaTeX
  \begin{align}
  
First Stage: Standard OLS \\
(1.1) X=Z\widehat{\gamma} \\ 
(1.2) Z'X=Z\widehat{\gamma} \\
(1.3) Z'X=Z'Z\widehat{\gamma} \\
(1.4) (Z'Z)^{-1}Z'X=(Z'Z)^{-1}Z'Z\widehat{\gamma} \\
(1.5) (Z'Z)^{-1}Z'X=\widehat{\gamma} \\


Second Stage, insert projection of Z on X ( \widehat{X} = Z \widehat{\gamma} ): \\
(2.1) Y= \widehat{X} \beta _{2SLS} \\
(2.2) \widehat{X}' Y= \widehat{X}' \widehat{X} \beta _{2SLS} \\
(2.3) ( \widehat{X}' \widehat{X})^{-1} \widehat{X}' Y=  ( \widehat{X}' \widehat{X})^{-1}  \widehat{X}' \widehat{X} \beta _{2SLS} \\
(2.4) ( \widehat{X}' \widehat{X})^{-1} \widehat{X}' Y=  \beta _{2SLS} \\
  
From (1.5) we insert  \widehat{X} = Z \widehat{\gamma} = Z (Z'Z)^{-1}Z'X \\
and  \widehat{X}' = \widehat{\gamma}'Z = X'Z (Z'Z)^{-1}Z'. \\
(2.5) ( X'Z (Z'Z)^{-1} Z'Z (Z'Z)^{-1} Z'X )^{-1}  X'Z (Z'Z)^{-1} Z'Y=  \beta _{2SLS} \\
(2.6) ( X'Z  (Z'Z)^{-1} Z'X )^{-1}  X'Z (Z'Z)^{-1} Z'Y=  \beta _{2SLS} \\

We will first derive the IV estimator: \\
(3.1) Y = X \beta _{IV} \\
(3.2) Z'Y = Z'X \beta _{IV} \\
(3.3) (Z'X)'Z'Y = (Z'X)'Z'X \beta _{IV} \\
(3.4) X'ZZ'Y = X'ZZ'X \beta _{IV} \\
(3.5) (X'ZZ'X)^{-1} X'ZZ'Y = (X'ZZ'X)^{-1} X'ZZ'X \beta _{IV} \\
(3.6) (X'ZZ'X)^{-1} X'ZZ'Y =  \beta _{IV} \\


To show that IV is same as the 2SLS we take an addition step after (3.2). Somewhat artificially we multiply both sides by X' Z (Z'Z)^{-1}  \\
(3.3b) X' Z (Z'Z)^{-1} Z'Y = X' Z (Z'Z)^{-1} Z'X \beta _{IV} \\
(3.4b) (X' Z (Z'Z)^{-1} Z'X)^{-1} X' Z (Z'Z)^{-1} Z'Y \\
                              = (X' Z (Z'Z)^{-1} Z'X)^{-1} X' Z (Z'Z)^{-1} Z'X \beta _{IV} \\
(3.5d) (X' Z (Z'Z)^{-1} Z'X)^{-1} X' Z (Z'Z)^{-1} Z'Y = \beta _{IV} \\
                     
  \end{align}
  Which must equal (3.6) because we have done only equal operations to both sides and (3.5d) is the same as (2.6).  Thus the IV estimator is the same as the 2SLS estimator. 
\end{document}