#include <oxstd.h>
#import <maximize>
#include <oxprob.h>


class ar_ml() {
	// variables
	decl y, pi;
	decl counter;	
	decl beta;
	decl std_err;
	decl qml_std_err;
	decl smpl;
	
	decl I2d;
	decl Iop;

	// functions
	ar2(const beta, const funct, const avScore, const amHessian);
	ar3(const beta, const funct, const avScore, const amHessian);
	outer_product();
	ar_ml(const data);
		
	estimate();
	print();
};


ar_ml::ar_ml(const y) {
	y = data;
	pi = acos(-1);
}


ar_ml::ar2(const beta, const funct, const avScore, const amHessian) {
	  // calculates exact log likelihood for an AR(1) process
	  // everything is vectorized, should run faster than the above version with a for loop
	  decl x, y1, u_sqr, logl, smpl;
	  x = (1 ~ lag0(y, 1))[1:][];
	 // print(x);
	  y1 = y[1:][];
	  smpl = sizer(y1);
	  logl = -0.5*((smpl+1)*log(2*pi) + smpl*log(beta[2][0]) + log(beta[2][0]/(1-beta[1][0]^2)));
	  u_sqr = ones(smpl, 1)'((y1 - x*beta[:1][]).^2);
	  logl = logl - 0.5*u_sqr/beta[2][0];
	  funct[0] = logl;	  
	  return 1;  
}


ar_ml::ar3(const beta, const funct, const avScore, const amHessian) {
	decl i = 1, logl = 0, error = 0;
	error = (y[counter][0] - beta[0][0]- beta[1][0] * y[counter - 1][0]);	
	funct[0] = -0.5*( log(2*pi) + log(beta[2][0]) + error^2 / beta[2][0]);
	return 1;
}


ar_ml::outer_product() {
	decl hh = zeros(3,3);
	decl cscore = zeros(3,1);
	decl smpl = sizer(y);
	for (counter = 1; counter < smpl; counter++) {
		Num1Derivative(ar3, beta, &cscore);
		hh = hh + cscore*cscore';
	}
	return hh;
}

ar_ml::estimate() {
	counter = 1;
	smpl = 50
	

main()
{
	counter = 1;
	decl smpl = 50;
	decl i, score, hess;
	decl ir, funct_value;
	decl e = rant(smpl + 50, 1, 20);
	//decl e = rann(smpl+50, 1);
	
	pi = acos(-1);
	y = zeros(smpl+50, 1);

	for(i = 1; i < smpl+50; i++) {                  // start this loop from 1 as have y(t-1) 
		y[i][0] = 1 + 0.3 * y[i-1][0] + e[i];
	}
	y = y[50:][];		
	decl beta = 0.1 + zeros(3, 1);
	decl at = ones(1,1);



	ir = MaxBFGS(ar2, &beta, &funct_value, 0 , TRUE);
	print("\nbeta: " , beta);

	Num2Derivative(ar2, beta, &hess);
  	print("std. errors as given by -E(inverse(information) " , diagonal(-1*invert(hess)).^0.5);
																		 

	decl iop = outter_product(beta); // information matrix as given by the outer product  
 	print("std. errors as given by outer product of derivatives " , diagonal(invert(iop)).^0.5 );
							 
	print("QMLE standard errors" , diagonal(invert(hess*invert(iop)*hess)).^0.5);
	
	savemat("c:/data/test.xls", y);
}