* File-Name:       OLS_in_MATA.do                                 
* Date:            January 23, 2009                                 
* Author:          Alex Herzog                                        
* Purpose:         This do-file gives a brief introduction to MATA
* Data Used:       tab9s (Jacobson and Dimock)
* Output File:     -
* Data Output:     -
* Machine:         Alex's MacBook                                     
                                                                    
clear                                                               
version 10.0                                                        
set more off
capture log close

* Loda dataset
use "tab9s.dta", clear

* OLS regression
regress cv92 lnce92

* Generate a constant (we will use this constant in MATA)
gen cons = 1

* ============================
* = Now let's deal with MATA =
* ============================
* Call MATA
mata

// Create the variables and take a look at them
// St-view is better than st_data
st_view(Y=., ., " cv92")
st_view(X=., ., ("cons","lnce92"))
Y
X

// First, find the OLS estimators
b = invsym(X'X)*X'Y
b

// Second, find the VC matrix
e  = Y - X*b
n  = rows(X)
k  = cols(X)
s2 = (e'e)/(n-k)
V = s2*invsym(X'X)
V

// This computes the standard errors of the estimators
se = sqrt(diagonal(V))

// Here's a nice presentation of results
(b, se)

// Now compute the test statistics
b:/se

// And p-values
2*normal(-abs(b:/se))

// Everything
(b, se, b:/se, 2*normal(-abs(b:/se)))


// Now, let's get some other stuff

// Create a scalar with the sum of squared residuals
// Check that it matches the results from the regression
RSS=e'e
RSS

// What are the RMSE? Just the RSS divided by N right?
RMSE=sqrt(RSS/n)
RMSE

// This is a cool matrix used to create deviations from means
st_view(i=., ., "cons")
iip=i*i'
mata describe iip
iip_over_n=iip/n
I309=I(309)

// The following matrix will help us to create deviations from means
Mo=I309-iip_over_n
mata describe Mo

// For instance, lets create the total sum of sum of squares
TSS=Y'*Mo*Y

// Which divided by N gives us the variance of Y
TSS/n

// Of course, we can also define a scalar and get means and variances of matrices
Y_mean=mean(Y)
Y_var=variance(Y)
Y_mean_var=meanvariance(Y)

// OK, having done this let's create R squared
// remember that R^2=1-RSS/TSS which we have
R2=1-RSS/TSS

// Now create the F-test: we will do loss of fit
st_view(Xr=., ., ("cons"))

// First, find the OLS estimators of restricted regression
br = invsym(Xr'Xr)*Xr'Y
br

// Second, create residuals from restricted regression
er  = Y - Xr*br

// Now the restricted RSS
RSSr=er'er
RSSr

// Get the elements of the F-test
J=1

// Compute it
F=((RSSr-RSS)/J)/(RSS/(n-k))

end
************************************
* now we are back in "normal Stata"
exit