/* GIBBS.G: sample program for ECON551 (Problem set 2)
**
**                             by Moto Shintani, February 1997
*/

new;
output file = gibbs.out reset;
format /rd 8,4;

/***(1)***/
y1 = 2;
y2 = 1;
rho = -.3;

n = 500;  /* sample size */
rndseed 234567;
z = rndn(n,2);  /* generate 2 independent standard normal r.v. */
theta1 = y1*ones(n,1) + z[.,1];  /* bivariate normal r.v. by Cholesky method*/
theta2 = y2*ones(n,1) + rho*z[.,1] +sqrt(1-rho^2)*z[.,2];

theta = theta1~theta2;
mean = meanc(theta);
vcov = vcx(theta);
output on;
print;
"Q4  MY NORMAL POSTERIOR";
print;
"           Sample mean     Sample covariance";
"******************************************************";
"(1) Direct method";
"<theta1>   ";;mean[1,1];;"       ";;vcov[1,1];
"<theta2>   ";;mean[2,1];;"       ";;vcov[2,1];;vcov[2,2];;
print;
output off;


/***(2)***/
t = 50;  /* # of iteration */

u = rndu(n,2);  /* generate 2 independent standard uniform r.v. */
theta1_0 = y1*ones(n,1) + (u[.,1] - .5); /* initial values */
theta2_0 = y2*ones(n,1) + (u[.,2] - .5);

i = 1;
do until i > n;
 z = rndn(t,2);  /* generate 2 independent standard normal r.v. */
 theta2_n = theta2_0[i]; /* theta1_0 is not used */
 j = 1;
 do until j > t;
  theta1_n = y1 + rho*(theta2_n - y2) + sqrt(1-rho^2)*z[j,1];
  theta2_n = y2 + rho*(theta1_n - y1) + sqrt(1-rho^2)*z[j,2];
  j = j + 1;
 endo;
 theta1[i] = theta1_n;
 theta2[i] = theta2_n;
 i = i + 1;
endo;

theta = theta1~theta2;
mean = meanc(theta);
vcov = vcx(theta);
output on;
print;
"(2) Gibbs sampler (500 initial values, 50 iterations)";
"<theta1>   ";;mean[1,1];;"       ";;vcov[1,1];
"<theta2>   ";;mean[2,1];;"       ";;vcov[2,1];;vcov[2,2];;
print;
output off;


/***(3)***/
t = 1000; /* # of iteration */

i = 1;
do until i > n;
 z = rndn(t,2);  /* generate 2 independent standard normal r.v. */
 theta2_n = theta2_0[i]; /* theta1_0 is not used */
 j = 1;
 do until j > t;
  theta1_n = y1 + rho*(theta2_n - y2) + sqrt(1-rho^2)*z[j,1];
  theta2_n = y2 + rho*(theta1_n - y1) + sqrt(1-rho^2)*z[j,2];
  j = j + 1;
 endo;
 theta1[i] = theta1_n;
 theta2[i] = theta2_n;
 i = i + 1;
endo;

theta = theta1~theta2;
mean = meanc(theta);
vcov = vcx(theta);
output on;
print;
"(3) Gibbs sampler (500 initial values, 1000 iterations)";
"<theta1>   ";;mean[1,1];;"       ";;vcov[1,1];
"<theta2>   ";;mean[2,1];;"       ";;vcov[2,1];;vcov[2,2];;
print;
output off;



/***(4)***/
t = 2500; /* # of iteration */
theta1 = zeros(t,1);
theta2 = zeros(t,1);

z = rndn(t,2);  /* generate 2 independent standard normal r.v. */
theta2_n = 0; /* initial value, theta1_0 is not used */
j = 1;
do until j > t;
 theta1_n = y1 + rho*(theta2_n - y2) + sqrt(1-rho^2)*z[j,1];
 theta2_n = y2 + rho*(theta1_n - y1) + sqrt(1-rho^2)*z[j,2];
 theta1[j] = theta1_n;
 theta2[j] = theta2_n;
j = j + 1;
endo;

theta1 = trimr(theta1,t-n,0);
theta2 = trimr(theta2,t-n,0);

theta = theta1~theta2;
mean = meanc(theta);
vcov = vcx(theta);
output on;
print;
"(4) Gibbs sampler (1 initial values, 2500 iterations)";
"<theta1>   ";;mean[1,1];;"       ";;vcov[1,1];
"<theta2>   ";;mean[2,1];;"       ";;vcov[2,1];;vcov[2,2];;
print;
output off;



end;