/** ECONOMETRICS 1 - SECOND EXAMPLE PROGRAM **/

/* simulation_example.txt */

/* in comparison to the first monte carlo exercise, we change the
   distributional assumption regarding the disturbance term. instead
   of a standard normal, the disturbance is a mixture of a normal and 
   a negative exponential
    random variable */

closeall;
output file = OLS_example.out reset;

/* y = 1 + 3x + epsilon */

/* set sample size */
n = 10;

/* x is uniform on 0,1 */
x = rndu(n,1);

/*generate the disturbance term*/

e1 = rndn(n,1);
e2 = -ln(1-rndu(n,1)) / .5 - 2;


/* weights applied to two random variables */
alpha = .2;
epsilon = alpha*e1 + (1-alpha)*e2;

bvec = 1 | 3;

xmat = ones(n,1) ~ x;

y = xmat * bvec + epsilon;

bhat = invpd(xmat'xmat)*xmat'y;

r = y - xmat*bhat;

s2 = r'r / (n - 2);

cov_bhat = s2*invpd(xmat'xmat);
se_bhat = sqrt(diag(cov_bhat));

/* compute t-test that b1 = 0 */
t = (bhat[2]-0)/se_bhat[2];
abst = abs(t);
pval = 2*cdftc(abst,n-2);



print "**************** OLS esimates *********************";
print "";
print "           bhat              se_bhat       ";
print (bhat ~ se_bhat);
print "t for H: b1 = 0" t;
print "p-value for this test stat" pval;


/* MONTE CARLO EXPERIMENT */

/* repeat on 100000 samples of size n - keeping x matrix fixed */

bmat = zeros(100000,2);
tstat = zeros(100000,1);

/* restrict DGP to be consistent with null hypothesis */
bvec_nullh = 1 | 0;

i = 1;
do until i gt 100000;
  new_epsilon = alpha*rndn(n,1) + (1-alpha)* (-ln(1-rndu(n,1)) / .5 - 2);
  new_y = xmat*bvec_nullh + new_epsilon;
  bhat_new = invpd(xmat'xmat)*xmat'new_y;
  bmat[i,1] = bhat_new[1];
  bmat[i,2] = bhat_new[2];
  newr = (new_y-xmat*bhat_new);
  s2 = newr'newr/(n-2);
  covmat = s2*invpd(xmat'xmat);
  tstat[i] = bhat_new[2]/(sqrt(covmat[2,2]));

  i = i+1;
endo;



print "**** means of b0 and b1 ******/";
print;
print meanc(bmat)';

/**
library pgraph;
graphset;
histp( tstat , 60 );
**/


/* compute p-value from the sampling distribution */
sim_pvalue = meanc( tstat .< -abs(t) .or tstat .> abs(t));
print "simulated pvalue" sim_pvalue;

/* if we don't know distribution of epsilon, can use residuals from unrestricted model
  and resample from those */

brep_mat = zeros(100000,2);
tstat2_vec = zeros(100000,1);

j = 1;
do until j gt 100000;
  sample_id = ceil(n*rndu(n,1));
  repr = r[sample_id,.];
  yrep = xmat*bvec_nullh + repr;
  brep = invpd(xmat'xmat)*xmat'yrep;
  brep_mat[j,.] = brep';
  newr = yrep - xmat*brep;
  s2 = newr'newr/(n-2);
  covmat = s2*invpd(xmat'xmat);
  tstat2 = brep[2] / sqrt(covmat[2,2]);
  tstat2_vec[j] = tstat2;
  
  j = j+1;
endo;


print "**** means of b0 and b1 ******/";
print;
print meanc(brep_mat)';


library pgraph;
graphset;
histp( tstat , 60 );


/* compute p-value from the sampling distribution */
sim_pvalue = meanc( tstat2_vec .< -abs(t) .or tstat2_vec .> abs(t));
print "simulated pvalue - bootstrap residuals" sim_pvalue;

/* bootstrapping 
   go back to original "sample",
    and redraw from that 

idea is to compute empirical distribution function of estimator or statistic*/

/* loop through this R times */
bigr = 20000;
brep_vec = zeros(bigr,2);

j = 1;
do until j gt bigr;
  sample_id = ceil(n*rndu(n,1));
  yrep = y[sample_id];
  xrep = xmat[sample_id,.];
  brep = invpd(xrep'xrep)*xrep'yrep;
  brep_vec[j,.] = brep';
  j = j+1;
endo;

histp(brep_vec[.,2],60);

/* compute probability that beta1 is 0 or less */

print "*** prob b1 is negative *****";
print meanc( brep_vec[.,2] .<= 0 );

end;