/******************************************************************
 A Mickey-Mouse example of the Gibbs sampler for a bivariate normal
with zero means and unit variances and correlation rho
******************************************************************/
rho = 0.99;      // correlation parameter
reps = 1000;     // replications in the multivariate simulation
s0 = 100;        // Burn-in for the Gibbs sampler
s1 = 1000;       // Actual draws used from the Gibbs sampler
/*******************************************************************
Method 1: Draw from a multivariate normal and induce the correlation
by appropriately normalizing the variables
*******************************************************************/
v = 1~rho|rho~1;    // Variance-Covariance matrix
p = chol(v);        // Cholesky decomposition s.t. p'p = v

z = rndn(reps,2);   // Generate the standard normals
y = z*p;            // Generate the correlation by appropriate rotation

"---------------------------------------------------";
"Variance-Covariance Matrix from Cholesky Simulation";
"---------------------------------------------------";
"Sample Mean: ";; meanc(y);
"Sample Variance-Covariance Matrix: ";
vcx(y);

/*******************************************************************
Use the Gibbs sampler to generate the multivariate normal draws
*******************************************************************/
y1 = zeros(s0+s1,1);        // Initialize the variables for the loop
y2 = zeros(s0+s1,1);

i = 2;
do while i le s0+s1;
    y2[i,1] = ((1 - rho^2)^0.5)*rndn(1,1) + rho*y1[i-1,1];
    y1[i,1] = ((1 - rho^2)^0.5)*rndn(1,1) + rho*y2[i,1];
    i = i+1;
endo;
ygs = y1~y2;
ygs = ygs[s0+1:s0+s1,.];
"---------------------------------------------------";
"Variance-Covariance Matrix from Gibbs Sampler";
"---------------------------------------------------";
"Sample Mean: ";; meanc(ygs);
"Sample Variance-Covariance Matrix: ";
vcx(ygs);

// Plot the histograms

if rows(y[.,1]) eq rows(ygs[.,1]);
begwind;
window(2,2,0);
setwind(1);
title("Histogram of Cholesky Monte Carlo for y1");
{ b,m,freq } = hist(y[.,1],20);
nextwind;
title("Histogram of Cholesky Monte Carlo for y2");
{ b,m,freq } = hist(y[.,2],20);
nextwind;
title("Histogram of Gibbs Sampler for y1");
{ b,m,freq } = hist(ygs[.,1],20);
nextwind;
title("Histogram of Gibbs Sampler for y2");
{ b,m,freq } = hist(ygs[.,2],20);
endwind;


endif;

/********************************************
Libraries
********************************************/
library pgraph;
graphset;