* Author: 	Le Wang
* Institute:  Population Center, University of Minnesota

		/*******************************************************************/
		/* Four Ways of Getting OLS Estimates in Mata			     */
		/*******************************************************************/

*-----------------------------------------------------------------------------------
* This program is to illustrate four different ways to calculate simple OLS 
* estimates in Mata. These cases seem quite trivial but give you some ideas 
* about how Mata works. Also, the time spent on calcluation for each case is 
* also reported so that you know in practice while different methods that give 
* you the same results may be quite different in terms of efficiency; in this 
* case, the difference is ignorable but may well be important in more complicated
* programs and larger samples.
*------------------------------------------------------------------------------------


webuse auto

mata
mata clear

	st_view(y=.,.,"weight",.)
	st_view(x=.,.,"length",.)

	/* Independent Variables	*/
	x=J(rows(x),1,1),x				/* Important trick to merge data		*/

	
	timer_clear()

/* Method 1			*/
timer_on(1)

		b=invsym(cross(x,x))*cross(x,y)

timer_off(1)


/* Method 2			*/
timer_on(2)

	b=invsym(x'*x)*(x'*y)

timer_off(2)


/* Method 3			*/
timer_on(3)

	void myols(todo, b, y, x, lnf, S, H)
	{
		lnf = (y - x*b') :* (y - x*b')	/*<--- Quiz: why?			*/

	}

	Sv = optimize_init()

	optimize_init_evaluator(Sv, &myols())

	optimize_init_evaluatortype(Sv, "v0")

	/* optimize_init_params(Sv, J(2,1,0))     <--- Quiz: why is it wrong?	*/

	optimize_init_params(Sv, J(1,2,0))
	optimize_init_which(Sv, "min")
	optimize_init_argument(Sv, 1, y)
	optimize_init_argument(Sv, 2, x)
	optimize(Sv)

timer_off(3)


/* Method 4			*/

timer_on(4)

	void myols1(todo, b, y, x, lnf, S, H)
	{
		lnf = (y - x*b')' (y - x*b')	/*<--- Quiz: Again, why?			*/

	}

	Sv = optimize_init()

	optimize_init_evaluator(Sv, &myols1())

	optimize_init_evaluatortype(Sv, "d0")    /*<---- Any difference here?		*/

	optimize_init_params(Sv, J(1,2,0))
	optimize_init_which(Sv, "min")
	optimize_init_argument(Sv, 1, y)
	optimize_init_argument(Sv, 2, x)
	optimize(Sv)

timer_off(4)

	
timer()	/* Report all the timers	*/
end


/* For comparison		*/
	reg weight length