%%HW 4 - Choosing the order for serial correlation test%%
clc
trial=5000;

n=1000;
e1=zeros(n,1); e2=zeros(n,1); e3=zeros(n,1);
test_1_reject=zeros(trial,3);
test_2_reject=zeros(trial,3);
beta1=zeros(trial,1); beta2=zeros(trial,1); beta3=zeros(trial,1);

for j=1:trial

for i=4:n
    e1(i,:)=randn(1,1);
    e2(i,:)=0.5*e2(i-1,:)+randn(1,1);
    e3(i,:)=0.4*e3(i-1,:)+0.3*e3(i-2,:)+0.2*e3(i-3,:)+randn(1,1);
end

e1=e1((n-99):n,:); e2=e2((n-99):n,:); e3=e3((n-99):n,:);
n=length(e1);

x=rand(n,1);
y1=0.5*x + e1; y2=0.5*x + e2; y3=0.5*x + e3;
beta1(j,:)=inv(x'*x)*x'*y1; 
beta2(j,:)=inv(x'*x)*x'*y2;  
beta3(j,:)=inv(x'*x)*x'*y3; 
resid1=y1-x*beta1(j,:); 
resid2=y2-x*beta2(j,:); 
resid3=y3-x*beta3(j,:);

%Test with 1 1ag

X1=[x(2:n,:) resid1(1:n-1,:)]; 
if (n-1)*(resid1(2:n,:)'*X1*inv(X1'*X1)*X1'*resid1(2:n,:))/(resid1(2:n,:)'*resid1(2:n,:))>3.84
    test_1_reject(j,1)=1;
end
X2=[x(2:n,:) resid2(1:n-1,:)]; 
if (n-1)*(resid2(2:n,:)'*X2*inv(X2'*X2)*X2'*resid2(2:n,:))/(resid2(2:n,:)'*resid2(2:n,:))>3.84
    test_1_reject(j,2)=1;
end
X3=[x(2:n,:) resid3(1:n-1,:)]; 
if (n-1)*(resid3(2:n,:)'*X3*inv(X3'*X3)*X3'*resid3(2:n,:))/(resid3(2:n,:)'*resid3(2:n,:))>3.84
    test_1_reject(j,3)=1;
end
%Test with 3 1ag
X1=[x(4:n,:) resid1(1:n-3,:) resid1(2:n-2,:) resid1(3:n-1,:)]; 
if (n-3)*(resid1(4:n,:)'*X1*inv(X1'*X1)*X1'*resid1(4:n,:))/(resid1(4:n,:)'*resid1(4:n,:))>7.81
    test_2_reject(j,1)=1;
end
X2=[x(4:n,:) resid2(1:n-3,:) resid2(2:n-2,:) resid2(3:n-1,:)]; 
if (n-3)*(resid2(4:n,:)'*X2*inv(X2'*X2)*X2'*resid2(4:n,:))/(resid2(4:n,:)'*resid2(4:n,:))>7.81
    test_2_reject(j,2)=1;
end
X3=[x(4:n,:) resid3(1:n-3,:) resid3(2:n-2,:) resid3(3:n-1,:)]; 
if (n-3)*(resid3(4:n,:)'*X3*inv(X3'*X3)*X3'*resid3(4:n,:))/(resid3(4:n,:)'*resid3(4:n,:))>7.81
    test_2_reject(j,3)=1;
end
end