#include <oxstd.h>
#import  <lib/probimhof>

//	Two examples of the function ProbImhof, included with probimhof.ox.

main()
{
	decl A, B, x, m, S, F, n, I, C, j, q; 
//	The first example checks the 95% quantile of the central F(5,3) distribution.
	A = (unit(5)~zeros(5,3))|zeros(3,8); A = A/5;
	B = zeros(5,8)|(zeros(3,5)~unit(3)); B = B/3;
	x = 9.01; m = zeros(8,1); S = unit(8);
	F = ProbImhof(x, A, B, m, S);
	println(" Prob [ F(5,3) <= 9.01 ] = ", F);

	q = QuanImhof(F, A, B, m, S);
	println(" Quan[ Prob [ F(5,3) <= 9.01 ] ] = ", q);

//	The second example gives the probability that the OLS estimator r
//	in an AR(1) model with unit root is greater than 1 (50 observations)
	n = 50;
	I = unit(n);
	C = zeros(n,n);
	for (j=1; j<n; ++j)
	{
		C = C + lag0(I,j);
	}
	A = C'(I+C); B = C'C; x = 1; m = zeros(n,1); S = unit(n);
	F = ProbImhof(x, A, B, m, S);
	println(" Prob [ r > 1 ] = ", (1-F));

	q = QuanImhof(F, A, B, m, S);
	println(" Quan[ Prob [ r > 1 ] ] = ", q);

	x = 0.9;
	F = ProbImhof(x, A, B, m, S);
	println(" Prob [ r > 0.9 ] = ", (1-F));

	q = QuanImhof(F, A, B, m, S);
	println(" Quan[ Prob [ r > 0.9 ] ] = ", q);
}
