new;
/* Sample GAUSS program for Econ 551b, PS1, Question 2 */
/* Hiu Man Chan, Yale University, January 1999 */

np=3; /* number of unknown parameters */
nsim=500; /* number of simulations */
n=1000; /* sample size in each simulation */
let beta[3,1]=1 0.5 1; /* true parameter value */

rndseed 234;

beta_ols=miss(0,0).*zeros(nsim,np); /* initialization to collect estimates */
i=1; do until i>nsim;
    ep=rndn(n,1);
    x0=ones(n,1);
    x1=rndn(n,1)-1;
    x2=0.5*rndn(n,1)+1;
    x=x0~x1~x2;
    y=exp(x*beta)+ep;
    beta_ols[i,.]=(inv(x'*x)*(x'*y))';
 
    i=i+1;
endo;

/* Write to output */
output file="ps1q2.out" reset;
"Econ 551b, PS1, Question 2";
"";"Results from Monte Carlo Experiment :";
"";"Summary Statistics of OLS estimates";
"Mean    Standard Deviation   Minimum  Maximum";    
meanc(beta_ols)~stdc(beta_ols)~minc(beta_ols)~maxc(beta_ols);
beta_max=ceil(maxc(beta_ols*20))/20;
beta_min=floor(minc(beta_ols*20))/20;
        
library pgraph;
_pdate="";
   
title("Distribution of beta_0");
ylabel("f(x)");xlabel("x");
   
br=seqa(beta_min[1],0.1,(beta_max[1]-beta_min[1])/0.1+1);

#IFUNIX
let v = 40 40 640 480 0 0 1 6 15 0 0 2 2;
wxy = WinOpenPQG(v,"XY Plot","XY");
call WinSetActive(wxy);
#ENDIF

{beta0_b,beta0_m,beta0_f}=histp(beta_ols[.,1],br);
   
#IFUNIX
call WinSetActive(1);
#ENDIF

title("Distribution of beta_0");
ylabel("f(x)");xlabel("x");

#IFUNIX
let v = 40 40 640 480 0 0 1 6 15 0 0 2 2;
wxy = WinOpenPQG(v,"XY Plot","XY");
call WinSetActive(wxy);
#ENDIF

xy(beta0_m,beta0_f);
   
#IFUNIX
call WinSetActive(1);
#ENDIF

title("Distribution of beta_1");
ylabel("f(x)");xlabel("x");
   
br=seqa(beta_min[2],0.1,(beta_max[2]-beta_min[2])/0.1+1);

#IFUNIX
let v = 40 40 640 480 0 0 1 6 15 0 0 2 2;
wxy = WinOpenPQG(v,"XY Plot","XY");
call WinSetActive(wxy);
#ENDIF

{beta1_b,beta1_m,beta1_f}=histp(beta_ols[.,2],br);
   
#IFUNIX
call WinSetActive(1);
#ENDIF

title("Distribution of beta_1");
ylabel("f(x)");xlabel("x");
   
#IFUNIX
let v = 40 40 640 480 0 0 1 6 15 0 0 2 2;
wxy = WinOpenPQG(v,"XY Plot","XY");
call WinSetActive(wxy);
#ENDIF

xy(beta1_m,beta1_f);
   
#IFUNIX
call WinSetActive(1);
#ENDIF

title("Distribution of beta_2");
ylabel("f(x)");xlabel("x");
   
br=seqa(beta_min[3],0.1,(beta_max[3]-beta_min[3])/0.1+1);

#IFUNIX
let v = 40 40 640 480 0 0 1 6 15 0 0 2 2;
wxy = WinOpenPQG(v,"XY Plot","XY");
call WinSetActive(wxy);
#ENDIF

{beta2_b,beta2_m,beta2_f}=histp(beta_ols[.,3],br);
   
#IFUNIX
call WinSetActive(1);
#ENDIF

title("Distribution of beta_2");
ylabel("f(x)");xlabel("x");
   
#IFUNIX
let v = 40 40 640 480 0 0 1 6 15 0 0 2 2;
wxy = WinOpenPQG(v,"XY Plot","XY");
call WinSetActive(wxy);
#ENDIF

xy(beta2_m,beta2_f);
   
#IFUNIX
call WinSetActive(1);
#ENDIF