/* A MONTE CARLO EXPERIMENT FOR THE GENERALIZED REGRESSION MODEL */

/* Create Data and Define Variables */
nobs = 20;
k = 1;
x = seqa(1,1,nobs);
rho = 0.9;
beta = 1;
sigma2 = 0.0625;

/* Create Variance-Covariance Matrix of Disturbances */
omega = eye(nobs);
for i (1,nobs,1);
	for j (1,nobs,1);
		omega[i,j] = rho^(abs(i-j))/(1-rho^2);
	endfor;
endfor;
invomega = inv(omega);

/* Variances of Estimators */
varbetahat = sigma2*inv(x'*invomega*x);
varb = sigma2*(inv(x'*x)*x'*omega*x*inv(x'*x));
exps2 = sigma2*sumc(diag(omega*(eye(nobs)-x*inv(x'*x)*x')))/(nobs-k);

/* Monte Carlo Experiment */

/* Storage Matrices */
loops = 5000;
betahat = zeros(loops,1);
b = zeros(loops,1);
s2star = zeros(loops,1);
s2 = zeros(loops,1);
vcvbetahat = zeros(loops,1);
vcvb = zeros(loops,1);

for i (1,loops,1);
		
	/* Create Artificial Data */
	v = zeros(nobs,1);
	u = rndn(nobs,1);
	v[1,1] = sqrt(sigma2/(1-rho^2))*u[1,1];
	for j (1,nobs-1,1);
		v[j+1,1] = rho*v[j,1] + sqrt(sigma2)*u[j+1,1];
	endfor;
	y = x + v;
	
	/* Create Statistics and Store in Matrices */
	betahat[i,1] = inv(x'*invomega*x)*x'*inv(omega)*y;
	b[i,1] = inv(x'*x)*x'*y;
	s2star[i,1] = (y - x*betahat[i,1])'*invomega*(y - x*betahat[i,1])/(nobs-k); 
	s2[i,1] = (y - x*b[i,1])'*(y - x*b[i,1])/(nobs - k);
	vcvbetahat[i,1] = s2star[i,1]*varbetahat/sigma2;
	vcvb[i,1] = s2[i,1]*inv(x'*x);
endfor;

/* Mean Results Across Simulations */
print "Average betahat = " meanc(betahat);
print;
print "Average b = " meanc(b);
print;
print "Actual sigma2 = " sigma2;
print;
print "Average s2star = " meanc(s2star);
print;
print "Expected value of s2 = " exps2;
print;
print "Average s2 = " meanc(s2);
print;
print "Actual var(betahat) = " varbetahat;
print;
print "Average Estimate of var(betahat) = " meanc(vcvbetahat);
print;
print "Actual var(b) = " varb;
print;
print "Average Estimate of var(b) = " meanc(vcvb);