p56 mylogit.ox

时间序列分析中著名的OxMetrics软件包

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

static decl y, x;

// We first need to generate data that follow a logit model
// with the supplied parameters. This Data Generator Process (dgp)
// is called the logistic function. Note how the binary 0/1 variable
// is generated so that it really follows the logit model.
// n is the number of observations.

logit_dgp(const n, const theta)
{
	decl k, x, y;
	
	k = sizer(theta); //sizer returns the number of rows of theta
	x = ones(n,1) ~ rann(n,k-1);

	y = (exp(x*theta) ./ (1 + exp(x*theta)) .> ranu(n,1));

// Note that the following dgp is equivalent 
// y = (1 ./ (1 + exp(-x*theta)) .> ranu(n,1));
	
	return y ~ x;
}

// We then have to define the log-likelihood function of this simple logit model.

logit(const theta, const obj, const score, const hess)
{
	decl prob, log_density, derivatives, ok;

	prob = 1 ./ (1 + exp(-x*theta));
	log_density = y .* log(prob) + (1 - y) .* log(1 - prob);

    obj[0] = meanc(log_density);  // average log-density

	return 1; // MaxBFGS needs the procedure to return 0 or 1
}

main()
{
  decl n, data, theta, theta_sim, error_code, obj_value, dfunc;
  	n = 1000; // starting values
	print("theta: ", theta_sim = 1|3);
	data = logit_dgp(n, theta_sim);
	println(data);

y=data[][0];
x=data[][1:];
obj_value = <0;0>;
MaxBFGS(logit, &obj_value, &dfunc, 0, TRUE);
print("function value is ", dfunc, 
        " at parameter value:", obj_value');

}