function [bhat,sigma,vcov1,vcov2,vcov1full]=varest(z,p,isintercept);
[bigt,n]=size(z);
y=z(p+1:bigt,:);
x=ones(bigt-p,1);
for j=1:p; x=[x z(p+1-j:bigt-j,:)]; end;
if isintercept==0; x=x(:,2:end); end;
bhat=y'*x*inv(x'*x);
res=y-(x*bhat');
sigma=(res'*res)/(bigt-(n*p)-1);
vcov1=kron(inv(x'*x),sigma);
d=duplication(n);
dplus=inv(d'*d)*d';
vcov2=2*dplus*kron(sigma,sigma)*(dplus');
vcov1full=vcov1;
vcov1=sqrt(diag(vcov1)); vcov2=sqrt(diag(vcov2));



