//////////////////////////////////////////
//Finite Lags
cls;
//n here is the number of observations
n = 100;
lag=5;
/* Data-Generating Process */
x1=ones(n,1);
x2 =rndu(n,1);
x3p=x2[1:n-1,1];
x3=rndu(1,1)~x3p';
x4p=x3[1,1:n-1];
x4=rndu(1,1)~x4p;
x5p=x4[1,1:n-1];
x5=rndu(1,1)~x5p;
x6p=x5[1,1:n-1];
x6=rndu(1,1)~x6p;
x7p=x6[1,1:n-1];
x7=rndu(1,1)~x7p;
x8p=x7[1,1:n-1];
x8=rndu(1,1)~x8p;
x9p=x8[1,1:n-1];
x9=rndu(1,1)~x9p;
x10p=x9[1,1:n-1];
x10=rndu(1,1)~x10p;
x11p=x10[1,1:n-1];
x11=rndu(1,1)~x11p;
x12p=x11[1,1:n-1];
x12=rndu(1,1)~x12p;
x13p=x12[1,1:n-1];
x13=rndu(1,1)~x13p;
x14p=x13[1,1:n-1];
x14=rndu(1,1)~x14p;
x15p=x14[1,1:n-1];
x15=rndu(1,1)~x15p;
x=x11'~x10'~x9'~x8'~x7'~x6';
//Coefficients
//these are our true values for data creation
beta = {6,5,4,3,2,1};
//creating data
y=x*beta+rndu(1,1)-.5;
//OLS estimates
b=inv(x'x)*(x'y);
print "Arithmetic Lag";
print "OLS estimate";
print b;
//Errors
eols=y-x*b;
ee=(eols'eols)/(100-6);
print"OLS S^2";
print ee;

//Estimating the Arithmetic Log using the procedure from
//Judge et al - basically an Almon lag with 1 degree.

//Creating our Z matrix, which is our data, 
//weighted by position. 
z = zeros(n,1);
zpart=zeros(6,1);
for j (1,n,1);
for i (0,5,1);
zpart[i+1,1]=(6-i).*x[j,i+1];
endfor;
z[j,1]=sumc(zpart);
endfor;

//Estimating alpha
//alpha is the slope of the weights on the lags
//it is constant because arithmetic lags are linear
//we do this by regressing y on z
alpha=(y'z)*inv(z'z);
print "Correct number of lags alpha, when the true is alpha=1";
print alpha;
//Backing out coefficient estimates using the definition of alpha
//where b=H*alpha
be=zeros(6,1);
for i (0,5,1);
//Estimate the betas using the equation betai=(n+1-i)alpha where (n+1-i)=H
be[i+1,1]=(6-i)*alpha;
endfor;
print "Correct lags coefficient estimate";
print be;
e=y-x[.,1:6]*be;
ee2=(e'e)/(100-2*5-1);
print"Correct number of lags squared errors";
print ee2;

//Arithmetic Lag where the number of lags is incorrect (guessing n=3 when the true is n=5)
//creating z matrix
z = zeros(n,1);
for j (1,n,1);
zpart=zeros(4,1);
for i (0,3,1);
zpart[i+1,1]=(4-i).*x[j,i+1];
endfor;
z[j,1]=sumc(zpart);
endfor;

//Estimate alpha 
// by regressing y on z
alpha=(y'z)*inv(z'z);
print "alpha where we incorrectly guess 3 lags";
print alpha;
//Coefficient estimates
be=zeros(4,1);
for i (0,3,1);
be[i+1,1]=(4-i)*alpha;
endfor;
print "estimate with too few lags";
print be;
e=y-x[.,1:4]*be;
ee2=(e'e)/(100-2*3-1);
print"S^2 when too few lags included";
print ee2;
print;
print;

/////////////////////////////////////////////
///////////Almon Polynomial//////////////////
//Almon lags require fitting both the number of lags and the degree
//of the polynomial, we can never have a higher polynomial degree
//than the number of lags. If the number of lags and degrees are
//known the process is simple
//degree
q=3;
//lags
n=5;
//storage
alpha=zeros(q+1,1);
H=zeros(n+1,q+3);
alpha={6,4,1,2};
print"Almon Lag";

//Filling our H matrix, which is the weight on alpha. By assumption beta=H*alpha
for i (1,n+1,1);
H[i,1]=1;
endfor;
for i (1,n+1,1);
H[i,2]=i-1;
endfor;
for i (1,n+1,1);
H[i,3]=(i-1)^2;
endfor;
for i (1,n+1,1);
H[i,4]=(i-1)^3;
endfor;
for i (1,n+1,1);
H[i,5]=(i-1)^4;
endfor;
for i (1,n+1,1);
H[i,6]=(i-1)^5;
endfor;
//Creating Data 
xx=x11'~x10'~x9'~x8'~x7'~x6';
xx2=x11'~x10'~x9'~x8'~x7'~x6'~x5'~x4'~x3';
y2=xx*H[.,1:4]*alpha+10*rndn(1,1);
//Actual beta
beta2=H[1:n+1,1:q+1]*alpha;
//ols
bb=inv(xx'xx)*(xx'y2);
e=y2-xx*bb;
olse=(e'e)/(100-n-q);


//Estimating alpha when we know the lag and polynomial degree
//We begin by deriving our Z matrix where z=xH
zz=xx[.,1:n+1]*H[1:n+1,1:q+1];

//Estimate alpha with a regression of y on z,
ahat=inv(zz'zz)*(zz'y2);
print "      actual alpha     estimated alpha";
print alpha~ahat;

//Estimating beta
//Solve for coefficient estimates using
//the relation between b and H and alpha
best2=H[1:n+1,1:q+1]*ahat;
e=y2-xx*best2;
este=(e'e)/(100-n-q);



print "S^2 when lags and degrees are known";
print este;
print "     actual beta       estimated beta    ols beta";
print beta2~best2~bb;


//Almon Polynomial Lag when too few lags
//Theory predicts that the estimates are biased
//q=degree
q=3;
//n=lags
n=4;
//Estimating alpha
zz=xx[.,1:n+1]*H[1:n+1,1:q+1];
ahat2=inv(zz'zz)*(zz'y2);
//Estimating beta
best3=H[1:n+1,1:q+1]*ahat2;
e=y2-xx[.,1:n+1]*best3;
ele=(e'e)/(100-n-q);
print "biased coefficients --- too few lags";
print best3;


//Almon Polynomial Lag when too few degrees
//theory predicts that the estimates will be biased
//q=degree
q=2;
//n=lags
n=5;
//Estimating alpha
zz=xx[.,1:n+1]*H[1:n+1,1:q+1];
ahat4=inv(zz'zz)*(zz'y2);
//Estimating beta
best4=H[1:n+1,1:q+1]*ahat4;
e=y2-xx[.,1:n+1]*best4;
ede=(e'e)/(100-n-q);
print "biased coefficients -- too few degrees";
print best4;

//Almon Polynomial Lag when too many degrees
//theory predicts that the estimates will be biased
//q=degree
q=4;
//n=lags
n=5;
//Estimating alpha
zz=xx[.,1:n+1]*H[1:n+1,1:q+1];
ahat5=inv(zz'zz)*(zz'y2);
//Estimating beta
best5=H[1:n+1,1:q+1]*ahat5;
e=y2-xx[.,1:n+1]*best5;
ede2=(e'e)/(100-n-q);
print "unbiased but inefficient coefficients --- too many degrees";
print best5;



print "ols S^2           estimated S^2  S^2 wrong lag    S^2 low degree      S^2 high degree";
print olse~este~ele~ede~ede2;
print;
print;

///////////////////////////////////////////////////////////////////////////////////
print"Fitting an unknown degree";
//We start by choosing the maximum possible poly degree that we think is possible
//This can be done using intuition or economic theory
deg=9;
n=10;
q=3;
//storage

H=zeros(n+1,deg+1);
alpha={6,4,1,2};

//fill H matrix
H2=zeros(n+1,deg+1);
for i (1,n+1,1);
H2[i,1]=1;
endfor;
for i (1,n+1,1);
H2[i,2]=i-1;
endfor;
for i (1,n+1,1);
H2[i,3]=(i-1)^2;
endfor;
for i (1,n+1,1);
H2[i,4]=(i-1)^3;
endfor;
for i (1,n+1,1);
H2[i,5]=(i-1)^4;
endfor;
for i (1,n+1,1);
H2[i,6]=(i-1)^5;
endfor;
for i (1,n+1,1);
H2[i,7]=(i-1)^6;
endfor;
for i (1,n+1,1);
H2[i,8]=(i-1)^7;
endfor;
for i (1,n+1,1);
H2[i,9]=(i-1)^8;
endfor;
for i (1,n+1,1);
H2[i,10]=(i-1)^9;
endfor;
//For illustration purposes, we increase the number of true lags to 10. 
xx2=x15'~x14'~x13'~x12'~x11'~x10'~x9'~x8'~x7'~x6'~x5';
y3=xx2*H2[.,1:q+1]*alpha+rndn(1,1);
btrue=H2[.,1:q+1]*alpha;

//storage
uc=zeros(deg+1,1);
s2=zeros(deg+1,1);
b3=zeros(n+1,deg+1);

//Process for estimating alpha and beta
for i (0,deg,1);
//R matrix, which is the restriction matrix
R=(H2[1:n+1,1:deg+1-i]*inv(H2[1:n+1,1:deg+1-i]'H2[1:n+1,1:deg+1-i])*H2[1:n+1,1:deg+1-i]');
//Estimating alpha by regressing y on z, changing the degree each time
z2=xx2[.,1:n+1]*H2[1:n+1,1:deg+1-i];
ahat=inv(z2'z2)*(z2'y3);
//Estimating beta using the assumed relationship between beta into h and alpha
be=H2[.,1:deg+1-i]*ahat;
//Errors
e2=y3-xx2[.,1:n+1]*be;
sig2=e2'e2/(100-n-i);
s2[i+1,1]=sig2;
//Performing our test, decreasing the polynomial degree each loop
u=(be'R*(R*inv(xx2[.,1:n+1]'xx2[.,1:n+1])*R')*R*be)/((100-deg-i)*sig2);
uc[i+1,1]=u;
b3[.,i+1]=be;
endfor;
print "crit=2.87 for 3 lags";
print "F-Stat";
print uc;
print "Estimates, changing the degree each time";
print "The first F-Stat corresponds with the last column below";
print "The first acceptable null hypothesis is the best fitting degree, and the corresponding estimates";
print "Are found by counting from the right";
//Our true beta values are in the first column, after that we have our estimates for the coefficients
//decreasing the number of polynomial degrees each time, coming from the right.
print "True betas        estimates, reducing the number of degrees each loop by 1";
print btrue~b3;
print;
print;

//////////////////////////////////////////////////////////////////////////////////////////
print "Fitting an unknown number of lags and unknown degree";
print "critical value is 2.2 at correct number";
//Following Pagano and Hartley 1975
//First we estimate lag length, then polynomial degree 
//We start with longest lag we are willing to consider, xN
xz=x11'~x10'~x9'~x8'~x7'~x6'~x5'~x4'~x3'~x2;
//We then use a gram-schmidt decomposition on xN
RR=qr(xz);
QQ=xz*inv(RR);
//Q is orthonormal and R is upper-triangular squared
//The first k columns of Q spand the same linear space as the first k columns of XN

//Use the Gram Schmidt decomposition to weight our data.
//We do this by rewriting the regression equation as y=Q*alpha+e using the relations X=QR and alpha=RB 
//We then regress y on Q to estimate alpha
//Q is orthogonal, or Q'Q=I so the regression is 
alphahat=QQ'y2;
//we find beta from the definition above
bph=inv(RR)*alphahat;
//we now have an unbiased estimate for alpha for the longest lag considered
//we now look for the last lag length where the restriction on our polynomial degrees fit

e2=y2-xz*bph;
ss2=e2'e2/(100-n);
test1=zeros(11,1);
//We repeatedly run an F-test to find for which polynomial degrees the model fits
for i (1,10,1);
test1[i,1]=alphahat[i,1]*alphahat[i,1]/((i+1)*ss2);
endfor;
//the first F-statistic that fails to reject holds all of the necessary lags. Because of this 
//the lag before this is the last lag that has to be included.
print"estimates    F-Stat";
print bph~test1[1:10,1];

//////////////////////////////////////////////////////////
//Infinite Lags
//Geometric Lag
//We begin by manipulating our infinite lag model
//y=alpha*sum(lambda*xt-1)+e
//by using the properties of infinite sums to derive the equation
//y=alpha*x+lambda*yt-1+e*

//Coefficients
//The true alpha and lambda coefficients are
coeff={1,.5};
y=zeros(100,1);
//creating data
y[1,1]=10;
Z=zeros(99,2);
for j (1,99,1);
	Z[j,.]=x2[j+1,1]~y[j,1];
	y[j+1,1]=Z[j,.]*coeff+rndn(1,1);
endfor;
//Koyck Process
//We begin by assuming that the original disturbance before our derivation, e was not autcorrelated
//even though our new error term e* is autocorrelated

//we run a simple ols regression on our transformed equation to get estimates of both alpha and lambda.
est=inv(Z'Z)*(Z'y[2:100,1]);
//however, we know that these estimates are biased, because plim(est)=(coeff)+E_xx^(-1)V
//where E=population covariace matrix of regressors
//and V=-lambda*sigma^2

f=y[1:99,1]'*y[2:100,1];


//Koyck respond to this by substititing their estimates into the new estimator 
//where we subtract the term "stuff" that has a probability limit equal to the bias in the OLS estimation
//we substitute our estimates from the OLS into the equation, and then solve for the actual values. 
//its a nonlinear system and involves a quadratic expression
//Where lambdae=the estimated value of lambda, and lambdat=the true value of lambda,
//stuff=(lambdat*Ee*^2)/(1+lambdat*lambdae)
a=x2[2:100,1]'*x2[2:100,1];
b=x2[2:100,1]'*y[1:99];
c=x2[2:100,1]'*y[1:99];
d=y[1:99,1]'*y[1:99,1];
e=x2[2:100,1]'*y[2:100,1];
//f=y[1:99,1]'*y[2:100,1]+stuff;

M=zeros(2,2);
M[1,1]=a;
M[2,1]=b;
M[1,2]=c;
M[2,2]=d;
N=zeros(2,1);
N[1,1]=e;
N[2,1]=f;
bb=inv(M)*N;
//Examples are taken from "The Theory and Practice of Econometrics" Judge Griffiths Hill Lee 1980