@ program = converge.g  ( EC670 ) @

    new;
    library pgraph;
    graphset;

    output file = converge.out reset;

/* Determine values for the simulation */

    iter = 1000;

/* Data Generating and calculating the sample mean:
    sample size = n, repeat ! for "iter" times      */

    b_keep = zeros(iter,2);
    v_keep  = zeros(iter,2);
    v1_keep  = zeros(iter,2);
    t_keep  = zeros(iter,2);
    t1_keep  = zeros(iter,2);

       a = 10;
       b = 5;
      sig = 10;
       b0 = a | b;

   "iter = " iter;
   "a  b  sig = " a  b  sig;

rhos = {0.0, 0.5, 1.0};
nobs = {10, 100, 1000, 3000};

rr=1;
do while rr <= rows(rhos);

  rho = rhos[rr];

  ""; "------------------"; 

nn=1;
do while nn <= rows(nobs);

  n = nobs[nn];

    if n > 1000; iter = 100; endif;

    i=1;
    do while i <= iter;

       e = rndn(n,1)*sig; 
       x = zeros(n,1); x[1] = rndn(1,1);
       y = zeros(n,1); y[1] = a + b * x[1] + e[1];
    
       j=2;    
        do while j <= n;
         x[j] = rho * x[j-1] + rndn(1,1);
         y[j] = a + b * x[j] + e[j];
         j = j + 1;
       endo;

        x = ones(n,1) ~ x;
        b1 = inv(x'x)*x'y;
        b_keep[i,.] = b1';
        r = y - x*b1;  
        sse = r'r/(n-2);
        var = sse * inv(x'x);
        v_keep[i,.] = diag(var)';

        if rho < 1;  v1_keep[i,.] = n * diag(var)';  endif;
        if rho == 1;  v1_keep[i,.] = n^2 * diag(var)';  endif;

        t_stat = (b1 - b0) ./ sqrt(diag(var));
        t1_stat = b1 ./ sqrt(diag(var));
        t_keep[i,.] = t_stat';
        t1_keep[i,.] = t1_stat';

        i = i + 1;
    endo;

    "";   
    format /rd/m1 6,1;
    " *** ";; "rho = ";; rho;; "   n = ";; n;     
    
    format /rd/m1 6,4;
    "est of b       = ";; meanc(b_keep)';; "   ";; 
    "var of b       = ";; meanc(v_keep)'; 

    "t_stat of b-b0 = ";; meanc(t_keep)';;  "   ";; 
    "t_stat of b    = ";; meanc(t1_keep)';

     "var of r(n)*b or n*b = ";; meanc(v1_keep)';
  
nn = nn+1;
endo;

rr = rr + 1;
endo;