/* computation class 2
   illustration of inversion problems
   solution of nonlinear equatio0ns  */

output file = comp_class_2.out reset;

/* when discussing algorithm efficiency and accuracy, we will need to calculate timing
   in GAUSS we can use the following builtin procedures for timing issues  */


x = hsec;  /*returns time since midnight in 1/100 of a second */

/* lets first look at some ill-conditioning problems, just in a linear estimation context */

/* imagine in the population the following linear relationship holds 
   
    y = 1 + .2x(1) - .3x(2) + e;

where e is N(0,1), and ( x(1) x(2) ) have a joint normal distribution with mean 0 and covariance
   matrix ( 2 -1, -1 1) 

To examine the performance of various algorithms for estimating this relationship using random samples,
we have to draw Monte Carlo samples from this population using pseudo-random number generators
we begin by generating sample draws for a sample of size N */

n = 1000000;

evec = rndn(n,1);

/* create bivariate normal matrix of regressors 

We can use the Cholesky decomposition to do this in a staightforward way,

for a MVN with mean vector mu and covariance matrix S, 

since S is PD can write  S = LL', where L is a lower triangular matrix

then 

x(1) = mu(1) + L(1,1)*rn(1)
x(2) = mu(2) + L(2,1)*rn(1) + L(2,2)*rn(2),

where x(i) is a standard normal rv draw */

s = {2 -1, -1 1};
lmat = chol(s)';

print lmat;

@ create random variables @

x = ones(n,3);
v = rndn(n,2);
x[.,2] = 0 + lmat[1,1] * v[.,1];
x[.,3] = 0 + lmat[2,1] * v[.,1] + lmat[2,2] * v[.,2];

print "sample covariance matrix" x'x/n ;

@ create y @

beta = { 1 , .2 , -.3 };

y = x*beta + evec;

/* compute regression coefficients */

/* time using inv operator */
h1 = hsec;
bhat1 = inv(x'x)*x'y;
h2 = hsec;
et1 = h2 - h1;

/* time using inv operator for pd matrices */
h1 = hsec;
bhat2 = invpd(x'x)*x'y;
h2 = hsec;
et2 = h2-h1;

/* other procedure */
h1 = hsec;
bhat3 = x'y / (x'x);
h2 = hsec;
et3 = h2-h1;

print (bhat1 ~ bhat2 ~ bhat3);
print;
print (et1 ~ et2 ~ et3);

/* now we investigate problems of ill-conditioning */
/* say that regression relationship is of the form

y = 1 + .00002 * x1 + 10000 * x2 + e

with x1 and x2 mean zero bivariate normal with covariance matrix (1000000 -30000, -30000 50) */

sigma = {1000000 -1000, -1000 5};

lmat = chol(sigma)';
print lmat;

beta_new = { 1 , .00002 , 10000 };

x2 = ones(n,3);
x2[.,2] = lmat[1,1] * v[.,1];
x2[.,3] = lmat[2,1] * v[.,1] + lmat[2,2] * v[.,2];

print "sample covariance matrix" x2'x2/n;

@ construct y @
y = x2 * beta_new + evec;

/* compute OLS estimates */

b2hat = invpd(x2'x2) * x2'y;

print "b2hat" b2hat;

end;