/* *********************************************************************************************************************************/
/* Estimates threshold  regression using least squares anmd two types of smoothed least squares modified 6/05/05                                                                                                */ 
/* *********************************************************************************************************************************/



new; 
cls;
screen on;
output file=thresholdrp.out reset ;		/** Where the rejection results are stored **/
screen off ;


t=round(100000*rndu(1,1)); 

cl = 0.001 ;
cu=1-cl;

load data=x;

n=rows(data);







y=data[2:n,1] ;
x=ones(n-1,1)~data[1:n-1,2] ;
n=n-1;


xt = x[.,2] ;
/* q=x[.,2] ;*/

/* q=seqa(1,1,n-1)/(n-1);*/
q=xt;
/* dat=y~x[.,2];*/





nb=200;

thetabb=zeros(nb,1) ;
thetasbb=zeros(nb,1) ;
thetaspbb=zeros(nb,1) ;
betbb=zeros(3,nb) ;
betsbb=zeros(3,nb) ;
betspbb=zeros(3,nb) ;
taus=zeros(nb,1) ;
tausp=zeros(nb,1) ;
sehb=zeros(3,1) ;
seshb=zeros(3,1) ;
sesphb=zeros(3,1) ;
tb=zeros(3,nb) ;
tsb=zeros(3,nb) ;
tspb=zeros(3,nb) ;




xname = "Constant "|"x"|"x1-" ;

k1=rows(x') ;
k2= rows(xt') ;
k=k1+k2;

n=rows(y) ;


/* np=n;*/
np=500;


par=zeros(np,1);

ms=zeros(k,k);



/*quantile(q,seqa(1/np,1/(np+1),np));*/
 
par = sortc(q,1);
s=stdc(q)*(log(n))*n^-0.5;


{ntheta,theta,thetas,thetasp,bet,bets,betsp,lik,liks,liksp,sbhomo,sbhetero,tse,tpse}=threshold(y,x,xt,q,par,s) ;


se= sbhomo[.,1] ;
ses= sbhomo[.,2] ;
sesp= sbhomo[.,3] ;


seh=sbhetero[.,1] ;
sesh=sbhetero[.,2] ;
sesph=sbhetero[.,3] ;


bs=bets[1:k1] ;
ds=bets[k1+1:k] ;
bsp=betsp[1:k1];
dsp=betsp[k1+1:k] ;
z = x~xt.*(q.<=theta) ;
zs= x~xt.*(ones(n,1)-ik((q-thetas*ones(n,1))./s)) ;
xtk=xt.*(ones(n,1)-ik((q-thetasp*ones(n,1))./s)) ;


resid= y - z*bet ;
resids= y - zs*bets ;
residsp= y - zs*betsp ;





/**************************************************************************************************************/
/*Outputs*/
/**************************************************************************************************************/

yname="returns" ;
xname = xname|"c" ;
qname = "Volatility" ;
est=bet|theta ;
se=seh|0.0000000000000000;

ff = "#*.*lG"~20~8;
"Unsmoothed Estimation";
print ("Dependent Variable:     " $+ yname);
print ("Regressor Variable:     " $+ qname);
print ("Threshold Variable:     " $+ qname);

  "Heteroskedasticity Correction Used";
"";
p = printfm("Variable"~"Estimate"~"St Error",0~0~0,ff);"";
"       __________________________________________";
p = printfm(xname~est~se,0~1~1,ff);"";
"";
"Observations                      " n;
"Degrees of Freedom                " (n-k);
"Sum of Squared Residuals          " meanc(resid.^2);
"";"";
"__________________________________________________";
"";"";


est=bets|thetas ;
se=sesh|tse;

ff = "#*.*lG"~20~8;
"Smoothed Estimation, Type I";
print ("Dependent Variable:     " $+ yname);
print ("Regressor Variable:     " $+ qname);
print ("Threshold Variable:     " $+ qname);

  "Heteroskedasticity Correction Used";
"";
p = printfm("Variable"~"Estimate"~"St Error",0~0~0,ff);"";
"       __________________________________________";
p = printfm(xname~est~se,0~1~1,ff);"";
"";
"Observations                      " n;
"Degrees of Freedom                " (n-k);
"Sum of Squared Residuals          " meanc(resids.^2);
"";"";
"__________________________________________________";
"";"";



est=betsp|thetasp ;
se=sesph|tpse;

ff = "#*.*lG"~20~8;
"Smoothed Estimation, Type II";

print ("Dependent Variable:     " $+ yname);
print ("Regressor Variable:     " $+ qname);
print ("Threshold Variable:     " $+ qname);

  "Heteroskedasticity Correction Used";
"";
p = printfm("Variable"~"Estimate"~"St Error",0~0~0,ff);"";
"       __________________________________________";
p = printfm(xname~est~se,0~1~1,ff);"";
"";
"Observations                      " n;
"Degrees of Freedom                " (n-k);
"Sum of Squared Residuals          " meanc(residsp.^2);
"";"";
"__________________________________________________";
"";"";


output off;
screen on;
closeall;
library pgraph ;

graphset ;

fonts("microb") ;
_pdate="" ;
_plegctl=0;





_paxes=1;
_pnum=2;
_plwidth=8;
/*
_paxht=0.5 ;
_pnumht=0.3;
_ptitlht=0.5;




_pltype={6 1 1 1 1};
_pstype=1 ;






*/

d=par~liks ;
d=sortc(d,1);
par=d[.,1] ;
liks=d[.,2] ;

_pcolor=4 ;

Title("Market Index Returns");

ylabel("Mean Sum of Squared Residuals");
xlabel("Threshold Parameter");
xy(par,liks);
/*xy(par./stdc(par),liks);*/
/* pause(5);*/
_pltype={6 6};
_plctrl={0 0 -1 } ;
ylabel("Returns");
xlabel("Volatility");



yfit=zs*bets ;
data=q~yfit~y;
data=sortc(data,1);

fx=data[1:ntheta-1,2]|miss(zeros(5,1),0) ;
fy=miss(zeros(ntheta+1,1),0)|data[ntheta+2:n,2] ;

xy(data[.,1],fx~fy~data[.,3]);



/**************************************************************************************************************/
/* The Main Procedure*/
/**************************************************************************************************************/


proc (14)=threshold(y,x,xt,q,par,s) ;

local ip,np,z,bet,resid,zs,bets,resids,lik,liks,liksp,xtk,m11,m12,m21,m22,ms,xtheta;
local betsp,bsp,dsp,bs,ds,residsp,residspsq,u,theta,us,thetas,usp,thetasp ;
local sbhomo, sbhetero,st1,st2, v,vsp,zsip,zsi,zi,vh,vsh,vsph,se,ses,sesp,seh,sesh,sesph;
local dresids,vv,hh,er,ee,ev,vs,hm,tse,tpse;
local k1,k2,k,n;
k1=rows(x') ;
k2= rows(xt') ;
k=k1+k2;
np=rows(par);
n=rows(y);
lik=zeros(np,1);
liks=zeros(np,1);
liksp=zeros(np,1);
ms=zeros(k,k);


for ip(1,np,1) ;

/**************************************************************************************************************/
/* Unsmoothed Estimator*/
/**************************************************************************************************************/


z = x~xt.*(q.<=par[ip]) ;

if det(z'z) .ne 0 ;
bet=inv(z'z)*z'y ;
resid= y - z*bet ;
else; 
resid= y;
endif;
lik[ip] = meanc(resid.^2) ; 

/**************************************************************************************************************/
/* Smoothed Estimator 1*/
/**************************************************************************************************************/


zs= x~xt.*(ones(n,1)-ik((q-par[ip]*ones(n,1))./s)) ;

if det(zs'zs) .ne 0 ;
bets=inv(zs'zs)*zs'y ;
resids= y - zs*bets ;
else;
resids=y;
endif;

liks[ip] = meanc(resids.^2) ; ;

/**************************************************************************************************************/
/*Smoothed Estimator 2*/
/**************************************************************************************************************/



xtk=xt.*(ones(n,1)-ik((q-par[ip]*ones(n,1))./s)) ;
m11 =x'x ;
m12 = x'xtk ;
m21=m12' ;
m22 = xt'xtk ;

ms[1:k1,1:k1]=m11;
ms[1:k1,k1+1:k]=m12;
ms[k1+1:k,1:k1]=m21;
ms[k1+1:k,k1+1:k]=m22;


if det(ms) .ne 0 ;
betsp=inv(ms)*(x'y|xtk'y);
bsp=betsp[1:k1];
dsp=betsp[k1+1:k] ;

residspsq= (y - x*bsp).^2  + ((xt*dsp) - 2*(y - x*bsp)).*(xtk*dsp);
else;
residspsq=y.^2;
endif;


liksp[ip] = meanc(residspsq) ; ;



/**************************************************************************************************************/
/* Return the Estimators*/
/**************************************************************************************************************/

endfor ;

u=lik~par~seqa(1,1,np);
u=sortc(u,1) ;
theta=u[1,2] ;
xtheta=u[1,3];


us=liks~par;
us=sortc(us,1) ;
thetas=us[1,2] ;


usp=liksp~par;
usp=sortc(usp,1) ;
thetasp=usp[1,2] ;



z = x~xt.*(q.<=theta) ;

bet=inv(z'z)*z'y ;


xtk=xt.*(ones(n,1)-ik((q-thetasp*ones(n,1))./s)) ;
m11 =x'x ;
m12 = x'xtk ;
m21=m12' ;
m22 = xt'xtk ;

ms[1:k1,1:k1]=m11;
ms[1:k1,k1+1:k]=m12;
ms[k1+1:k,1:k1]=m21;
ms[k1+1:k,k1+1:k]=m22;

betsp=inv(ms)*(x'y|xtk'y);

zs= x~xt.*(ones(n,1)-ik((q-thetas*ones(n,1))./s)) ;
bets=inv(zs'zs)*zs'y ;


/********************************

Standard Errors for Slopes 

**********************************/


bs=bets[1:k1] ;
ds=bets[k1+1:k] ;
bsp=betsp[1:k1];
dsp=betsp[k1+1:k] ;

resid= y - z*bet ;
resids= y - zs*bets ;
residsp= y - zs*betsp ;



zi=z.*(resid);
vh=zi'zi ;
zsi = zs.*(resids);
vsh = zsi'zsi;
zsip=(x~xtk).*(residsp);
vsph = zsip'zsip ;

vh=inv(z'z)*vh*inv(z'z);
vsh=inv(zs'zs)*vsh*inv(zs'zs);
vsph=inv(ms)*vsph*inv(ms);



v = meanc(resid.^2)*inv(z'z) ;
vs = meanc(resids.^2)*inv(zs'zs) ;
vsp =   meanc(residsp.^2)*inv(ms)*(zs'zs)*inv(ms) ;

se=sqrt(diag(v)) ;
ses=sqrt(diag(vs)) ;
sesp=sqrt(diag(vsp)) ;

seh=sqrt(diag(vh)) ;
sesh=sqrt(diag(vsh)) ;
sesph=sqrt(diag(vsph)) ;


/********************************

Standard Errors for Threshold Parameter

**********************************/



/* Smoothed Estimator 1 */

dresids = (xt*ds).*dk((q-thetas*ones(n,1))./s)./s ;

VV = sumc((resids.^2).*(dresids.^2)) ;
HH=sumc((dresids.^2)) ;

tse = sqrt(VV/HH^2) ;



/* Smoothed Estimator 2 */

er = y - x*bsp ;

ee = (  (xt*dsp).^2 - (2*xt*dsp).*er  ) ;


ev = ee.*dk((q-thetasp*ones(n,1))./s)./s ;
Vs =sumc(ev.^2) ;
Hm = sumc(ee.*ddk((q-thetasp*ones(n,1))./s))./(s.^2) ; 

tpse = sqrt(Vs/(Hm^2)) ;



sbhomo=se~ses~sesp;
sbhetero=seh~sesh~sesph;
st1=tse ;
st2=tpse ;



retp(xtheta,theta,thetas,thetasp,bet,bets,betsp,lik,liks,liksp,sbhomo,sbhetero,st1,st2);
endp;








/* This procedure computes the vector of univariate density estimates at evaluation points xe using data x*/

proc matkerd(xe,x,h) ;
local m,n,nne,K, K1,k2,k3,hl ;

n=rows(x) ;
nne=rows(xe) ;

hl=h ;
m = (xe*ones(1,n) - ones(nne,1)*x')./hl ; 
K1 = pdfn(m)./hl ;

K=(meanc(K1'))' ;

retp(K) ;
endp ;


/* standard bootstrap algorithm (vm) */

proc boot_vm(x);

	local n,y;

	n = rows(x);

	y  =  x[trunc(rndu(n,1)*n)+1,1];
y;
	retp(y);

endp;


/* standard pair bootstrap algorithm (vm) */

proc bootp_vm(x);

	local n,y;

	n = rows(x);

	y  =  x[trunc(rndu(n,1)*n)+1,.];



	retp(y);

endp;





fn ik(x) = cdfn(x) + x.*pdfn(x);
fn dk(x)= (2 - x.^2).*pdfn(x) ;
fn ddk(x) = (x.^3 -3).*pdfn(x) ;