function [muhat,Kmathat,sigmahat,tstat,tstatwhite,qvals,pvals,qvalswhite,pvalswhite]=ptv(Y,X,thetanull,tau);

% This function carries out nearly exactly inference in the following
% model:
% y(t) = mu1 + theta*x(t-1) + e1(t)
% x(t) = mu2 + rho*x(t-1) + e2(t)
% We assume the errors e1 and e2 are mean zero normals,
% possibly contemporaneously correlated, and independent
% across time.

muhat=zeros(2,1); Kmathat=zeros(2,2); sigmahat=zeros(2,2);
%%%%%%%%%%%%%%%%%%%%%%%%
% Construct point estimates
n=length(X);
% The y(t) on x(t-1) equation:
ymat=Y(2:n);
xmat=[ones(n-1,1),X(1:n-1)];
bols=inv(xmat'*xmat)*(xmat'*ymat);
muhat(1)=bols(1);
Kmathat(1,2)=bols(2);
r1=ymat-xmat*bols;
% The x(t) on x(t-1) equation:
ymat=X(2:n);
xmat=[ones(n-1,1),X(1:n-1)];
bols=inv(xmat'*xmat)*(xmat'*ymat);
muhat(2)=bols(1);
Kmathat(2,2)=bols(2);
r2=ymat-xmat*bols;

%%%%%%%%%%%%%%%%%%%%%%%%%
% estimate error covariance matrix
s1mle=sqrt(r1'*r1/(n-2));
s2mle=sqrt(r2'*r2/(n-2));
gammamle=(r1'*r2/(n-2))/(s1mle*s2mle);
sigmahat=[s1mle*s1mle,s1mle*s2mle*gammamle;s1mle*s2mle*gammamle,s2mle*s2mle];

%%%%%%%%%%%%%%%%%%%%%%%%
% We are going to test a series of nulls
nnull=length(thetanull);
tstat=zeros(nnull,1);
qvals=zeros(nnull,1);
pvals=zeros(nnull,1);
tstatwhite=zeros(nnull,1);
qvalswhite=zeros(nnull,1);
pvalswhite=zeros(nnull,1);
for inull=1:nnull;

	Ynull=Y(2:n)-thetanull(inull)*X(1:n-1);
	%%%%%%%%%%%%%
	% What is the t statistic?
	ymat=Ynull;
	xmat=[ones(n-1,1),X(1:n-1)];
	bols=inv(xmat'*xmat)*(xmat'*ymat);
	r1=ymat-xmat*bols;
	s1=sqrt(r1'*r1/(n-2));
	Xdemean=X(1:n-1)-mean(X(1:n-1));
	V=s1*s1/sum(Xdemean.^2);
	tstat(inull)=bols(2)/sqrt(V);

	% what about the heteroskedasticity-robust tstat?
	Xdemean=X(1:n-1)-mean(X(1:n-1));
	Vwhite=sum(Xdemean.*Xdemean.*r1.*r1)/(sum(Xdemean.*Xdemean))^2;
	tstatwhite(inull)=bols(2)/sqrt(Vwhite);

	% Constrained mle for rho
	s1=s1mle; s2=s2mle; gamma=gammamle;
	ymat=X(2:n)/s2-gamma*(Ynull)/s1;
	xmat=[ones(n-1,1),X(1:n-1)/s2];
	bols=inv(xmat'*xmat)*(xmat'*ymat);
	Khat=n*(bols(2)-1);
	% Third, hessian term
	Hstat=inv(xmat'*xmat);
	Hhat=(1.0/Hstat(2,2))/(n*n);
	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	% The sufficient statistics (plus info about the sample size) are
	SS=[exp(Khat/50);exp(-Hhat);abs(gamma);exp(-n/100)];

	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	% neural net pvals and quantiles for homoskedastic tstat
	psi =[    1.8702    3.7744   49.9034   -1.6534
    2.1040   -2.5565    2.9268   -1.0395  
   -3.4355   -1.9475  -52.7576    2.8437  
    0.5738   -0.8120   -5.0194   -0.2264  
    0.0119   -0.0262    4.4890    0.0084];

	Pmean =[   -1.7383
   	4.5693
   	2.7826
    -0.0007
   	3.9894	];

	Pstd =  [ 0.1746
   -0.4631
   -0.1955
    0.0210
   -0.4641];

	Xnn=[1;SS];
	G=ones(5,1);
	for j=2:5;
		G(j)=tanh(psi(:,j-1)'*Xnn);
	end;

	qvals(inull) = sign(gamma)*Pmean'*G + exp(Pstd'*G)*norminv(tau);
	pvals(inull) = 1 - normcdf((tstat(inull) - sign(gamma)*Pmean'*G)/exp(Pstd'*G));

	%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%%
	% the heteroskedasticity-robust statistic is asymptotically
	% equivalent to the regular t-statistic, under both stationary
	% asymptotics and the local-to-unity asymptotics. So for
	% simplicity I'll just use the same neural net for the
	% heteroskedasticity-robust statistic
	Xnn=[1;SS];
	G=ones(5,1);
	for j=2:5;
		G(j)=tanh(psi(:,j-1)'*Xnn);
	end;

	qvalswhite(inull) = sign(gamma)*Pmean'*G + exp(Pstd'*G)*norminv(tau);
	pvalswhite(inull) = 1 - normcdf((tstatwhite(inull) - sign(gamma)*Pmean'*G)/exp(Pstd'*G));
	
end;