clear;
randn('state',sum(100*clock));
rand('state',sum(100*clock));

% simulate data from the model
% y(t) = theta*x(t-1) + e1(t)
% x(t) = rho*x(t-1) + e2(t)
rho=1;
gamma=-.9;
theta=0;
K=[0,theta;0,rho];
SIGMA=[1,gamma;gamma,1];
R=(chol(SIGMA))';

% number of observations
n=12*50;

fprintf('\nSimulating %d monthly observations.',n);
fprintf('\nTrue parameters: theta is %g, gamma is %g, rho is %g.',theta,gamma,rho);

% simulate the data
DATA=zeros(2,n);
Z=randn(2,n);
E=R*Z;
DATA(1,:)=E(1,:);
for j=2:n;
	DATA(:,j) = K*DATA(:,j-1) + E(:,j);
end;    
Y=DATA(1,:)';
X=DATA(2,:)';

% test
thetanull=(-.05:.005:.05)';
taulow=.025;
[muhat,Kmathat,sigmahat,tstat,tstatwhite,qvals05,pvals,qvalswhite05,pvalswhite]=ptv(Y,X,thetanull,taulow);
tauhigh=.975;
[muhat,Kmathat,sigmahat,tstat,tstatwhite,qvals95,pvals,qvalswhite95,pvalswhite]=ptv(Y,X,thetanull,tauhigh);

fprintf('\nEstimated model:');
fprintf('\n y(t) = %g + %g * x(t-1) + e1(t)',muhat(1),Kmathat(1,2));
fprintf('\n x(t) = %g + %g * x(t-1) + e2(t)',muhat(2),Kmathat(2,2));
fprintf('\n Stddev(e1) = %g, Stddev(e2) = %g, Corr(e1,e2) = %g,',sqrt(sigmahat(1,1)),sqrt(sigmahat(2,2)),sigmahat(1,2)/sqrt(sigmahat(1,1)*sigmahat(2,2)));

fprintf('\nTesting nulls with homoskedastic tstats');
fprintf('\nNull\ttstat\t%g quantile\t%g quantile',taulow,tauhigh);
nnull=length(thetanull);
for inull=1:nnull;
	fprintf('\n%5.3f\t%5.3f\t%5.3f\t\t\t%5.3f',thetanull(inull),tstat(inull),qvals05(inull),qvals95(inull));
end;

fprintf('\nTesting nulls with heteroskedastic tstats');
fprintf('\nNull\ttstat\t%g quantile\t%g quantile',taulow,tauhigh);
nnull=length(thetanull);
for inull=1:nnull;
	fprintf('\n%5.3f\t%5.3f\t%5.3f\t\t\t%5.3f',thetanull(inull),tstatwhite(inull),qvalswhite05(inull),qvalswhite95(inull));
end;