** This procedure performs Tsay's (1989)
** TAR specification test
** and was written by Ming Chien Lo on July 13, 1999.
**
** Last updated: August 28, 1999.
**
** Input: z T x 1 data series
**
** p order of the autoregressive process
** d lag index for threshold variable Y{t-d}
**
** Output: F resulting F-statistics
** b starting observation index for the recursive regression
** h starting index of the arranged data set
**
** (both b & h are for calculation of d.f. in F-test)
*/
proc(3)=tsayuni(z,p,d);
local h, t, oz, ytd, y, x, i, y_, x_, Pm, bm, ym1, xm, xm1, ahat, beta, Dm1, Dm1t, Fd_, Fn_, F,
Km1, Pm1, bm1, b, ehat, ephat, omega, invXX, gs, zx, xtd, x1, yt, xt;
h=maxc(1|(p+1-d)); @ starting index of the arranged data set @
t=rows(z); @ total no. of observations @
ytd=z[h:t-d]; @ potential threshold @
yt=z[h+d:t]; @ truncated y series @
zx=shiftr(z',seqa(1,1,p),0); @ create independent variable matrix from raw data @
zx=zeros(p,p)|trimr(zx',p,0);
xt=zx[h+d:t,.]; @ trucated x series @
gs=sortind(ytd);@ index for sorted y{t-d} @
y=yt[gs]; @ arranged Y according to sorted y{t-d} @
x=xt[gs,.]; @ arranged X according to sorted y{t-d} @
b=floor(0.1*t)+p;@ set starting observation index for the recursive regression @
x=ones(rows(x),1)~x;
@ derive first beta_m for the first m cases according to b @
beta=zeros(t-d-h-b+1,p+1);
Dm1t=zeros(t-d-h-b+1,1);
ahat=zeros(t-d-h-b+1,1);
y_=y[1:b];
x_=x[1:b,.];
Pm=invpd(moment(x_,0));
bm=Pm*x_'y_;
@ recursive regression @
i=1;
do until i>t-d-h-b+1;
ym1=y[b+i];
xm1=x[b+i,.]';
Dm1=1+xm1'*Pm*xm1;
Km1=Pm*xm1/Dm1;
bm1=bm+Km1*(ym1-xm1'*bm);
Pm1=(eye(p+1)-Pm*(xm1*xm1')/Dm1)*Pm;

beta[i,.]=bm1';
Dm1t=Dm1;
ahat=ym1-xm1'*bm;
bm=bm1;
Pm=Pm1;
i=i+1;
endo;
@ calculate F-statistics @
ehat=ahat./sqrt(Dm1t);
x1=x[b+1:rows(x),.];
invXX=invpd(moment(x1,0));
omega=invXX*x1'ehat;
ephat=ehat-x1*omega;
Fn_=(sumc(ehat.*ehat)-sumc(ephat.*ephat))/(p+1);
Fd_=sumc(ephat.*ephat)/(t-d-b-p-h);
F=Fn_/Fd_;
retp(F,b,h);
endp;