* GIBBS.DO

log using gibbs.txt, text replace

* 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
  v = (1, $rho \ $rho, 1)  // Variance-Covariance matrix
  p = cholesky(v)          // Cholesky decomposition s.t. p'p = v
  p
  z = invnormal(uniform($reps,2))  // Generate the standard normals
  y = z*p'             // Generate correlation draws by appropriate rotation
  rows(y)
  mean(y,1)
  variance(y,1)
end

********** METHOD 2: GIBBS SANPLER **********

mata
  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)
  mean(ygs,1)
  variance(ygs,1)
  last20ygs = ygs[|($s1-19),1 \ ($s1),.|]
  last20ygs
end

********* CLOSE OUTPUT **********
log close
clear 
exit