---------------------------------------------------------------------------------
       log:  c:\TEACHING\bayeskoop\gibbs.txt
  log type:  text
 opened on:  27 Oct 2005, 17:18:18

. 
. * Program gibbs.do October 2005 for Stata version 9
. * Based on Oscar Jorda's Gauss program October 21 2005 
. 
. ********** OVERVIEW OF ASS1F05.DO **********
. * This program does 
. *   (1) Gibbs sampler for bivariate normal [0, I] 
. * For     rho = .5 S0 = S1 = 1000 okay
. * but for rho = .9 need S0 = S1 = 10000
. 
. ********** SETUP **********
. * set mem 5m
. set more off

. version 9

. 
. ********** GLOBALS USED BELOW **********
. 
. set seed 1234567       // So get same draws in future runs  

. global rho = 0.90      // Correlation parameter

. global reps = 10000    // Replications in the Cholesky draws

. global s0 = 10000      // Burn-in for the Gibbs sampler

. global s1 = $reps      // Actual draws used from the Gibbs sampler

. 
. ********** METHOD 1: USUAL WITH CHOLESKY DECOMPOSITION **********
. 
. mata
------------------------------------------------- mata (type end to exit) -----
:   v = (1, $rho \ $rho, 1)  // Variance-Covariance matrix

:   p = cholesky(v)          // Cholesky decomposition s.t. p'p = v

:   p
                 1             2
    +-----------------------------+
  1 |            1             0  |
  2 |           .9   .4358898944  |
    +-----------------------------+

:   z = invnormal(uniform($reps,2))  // Generate the standard normals

:   y = z*p'             // Generate correlation draws by appropriate rotation

:   rows(y)
  10000

:   mean(y,1)
                 1             2
    +-----------------------------+
  1 |  .0173902148   .0115941983  |
    +-----------------------------+

:   variance(y,1)
[symmetric]
                 1             2
    +-----------------------------+
  1 |  1.000220523                |
  2 |  .9016835667    1.00100506  |
    +-----------------------------+

: end
-------------------------------------------------------------------------------

. 
. ********** METHOD 2: GIBBS SANPLER **********
. 
. mata
------------------------------------------------- mata (type end to exit) -----
:   y1 = J($s0+$s1,1,0)  // Initialize the variables for the loop

:   y2 = J($s0+$s1,1,0)

:   i = 2

:   do {
>       y1[i,1] = ((1 - $rho^2)^0.5)*invnormal(uniform(1,1)) + $rho*y2[i-1,1]
>       y2[i,1] = ((1 - $rho^2)^0.5)*invnormal(uniform(1,1)) + $rho*y1[i,1]
>       i = i+1
>   } while (i <= $s0+$s1)

:   ygs = y1,y2

:   ygs = ygs[|($s0+1),1 \ ($s0+$s1),.|]  //Drop the burn-ins

:   rows(ygs)
  10000

:   mean(ygs,1)
                  1              2
    +-------------------------------+
  1 |   -.033560495   -.0301804643  |
    +-------------------------------+

:   variance(ygs,1)
[symmetric]
                 1             2
    +-----------------------------+
  1 |  1.028309198                |
  2 |  .9286431487   1.024096204  |
    +-----------------------------+

:   last20ygs = ygs[|($s1-19),1 \ ($s1),.|]

:   last20ygs
                   1              2
     +-------------------------------+
   1 |  -.3089312273   -.1268435449  |
   2 |   .6864120381     1.44286724  |
   3 |   1.607474633    1.335187689  |
   4 |   .9520909338    .6059896901  |
   5 |   .7316986634    .1628200431  |
   6 |   .3462271484    .3091596226  |
   7 |   .6720405225    -.134106668  |
   8 |   .0126563281   -.3471433691  |
   9 |  -.2305736159    .5314264302  |
  10 |   1.321107413    1.892644655  |
  11 |   1.910559619    1.776674083  |
  12 |   1.753229729    1.473484238  |
  13 |   1.995493843    2.020478429  |
  14 |   2.215578025    2.062270339  |
  15 |   1.916963858    1.869871988  |
  16 |   1.834126399    .4721176697  |
  17 |   .9345823275    .9865429189  |
  18 |   1.161096657    .9385897539  |
  19 |   .2768269394    .8405649927  |
  20 |   1.254689905    1.804056879  |
     +-------------------------------+

: end
-------------------------------------------------------------------------------

. 
. ********* CLOSE OUTPUT **********
. log close
       log:  c:\TEACHING\bayeskoop\gibbs.txt
  log type:  text
 closed on:  27 Oct 2005, 17:18:18
-------------------------------------------------------------------------------