new;
cls;
library pgraph;


/* This program computes the order of summability estimator and its
corresponding subsampling confidence intervals. The partial demeaning 
case that considers the presence of a constant term is considered. */

/* Load or Generate Data */

load fedfunds[682,1]=C:\Users\Vanessa\Dropbox\Vanessa\GaussV\PhD\data\taylorrule\fedfunds.txt;

z0t=fedfunds;
T=rows(z0t);
b=int(sqrt(T))+1;

/* Deterministics */

/* Constant Case: Partial Demeaning */

zt=zeros(T,1);
van=1;
do until van>T;
zt[van]=z0t[van]-(sumc(z0t[1:van])/van);
van=van+1;
endo;
zt=zt[2:T];

T=rows(zt);
b=int(sqrt(T))+ 1; 


/* Estimator of the Full Sample */

xdtn=ln(seqa(1,1,T));
ydtn=ln(cumsumc(zt).^2);
ydtn=ydtn-ydtn[1];

bmcon=inv(xdtn'*xdtn)*(xdtn'*ydtn);

"Full Sample Estimator";
"beta=";bmcon;
"delta=";(bmcon-1)/2; ?;


/******* Subsampling samples generator *******/

numsubsam=T-b+1;

subsamples=zeros(b,T-b+1);

i=1;
do while i <= T-b+1;

subsamples[.,i]=zt[i:b+i-1];

i=i+1;
endo;

/*** Subsamples Estimation ***/


Nn=cols(subsamples);
bmco=zeros(Nn,1);

j=1;
do until j>Nn;

xt=ln(seqa(1,1,b));
yt=ln(cumsumc(subsamples[.,j]).^2);
yt=yt-yt[1];

bmcosub=inv(xt'*xt)*(xt'*yt);

bmco[j]=ln(b)*(abs(bmcosub-bmcon));

j=j+1;
endo;

/*** Empirical CDF ***/

{y,c}=pempiric(bmco);

/*** Critical Values ***/

alphasig=0.05;
confidence=1-alphasig;
confivector=ones(rows(c),1).*confidence;
difconfic=abs(c-confivector);
position=minindc(difconfic);
cnb1malphahalph=y[position,1];

/*** Confidence Intervals ***/

Ilow=bmcon-((1/ln(T))*cnb1malphahalph);
Iup=bmcon+((1/ln(T))*cnb1malphahalph);
"Subsampling Interval";
"betas=";
Ilow~Iup;
"deltas=";
((Ilow-1)/2)~((Iup-1)/2);


/*
**  empiric.src
**
**  Version 1.0  16/05/2000 (C) Copyright 2000 by Rainer Schlittgen
**  All rights Reserved.
**
**   Frequencies, empirical cumulative distribution function
**   Use the GAUSS procedure QUANTILE for the inverse empirical cdf
**
**  Format
**  ----------------------------------------------------------------------
**  {y,f} = dempiric(x)                                         
**  {y,c} = pempiric(x)                                         
**
**  Input
**  ----------------------------------------------------------------------
**  x...(k,m)-matrix of observed (univariate) values
**
**  Output
**  -------------------------------------------------------------------
**  y...(n,1)-vector, of ordered distinct values in x
**  f...(n,1)-vector, of rel. frequencies at y
**  c...(n,1)-vector, of cumulated rel. frequencies at y
**
**********************************************************************/

proc (2) = dempiric(x);
 local f,y;
  x = VEC(x);
  y = SORTC(UNIQUE(x,1),1);
  f = COUNTS(x,y)/ROWS(x);
 retp(y,f);
endp;

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

proc (2) = pempiric(x);
 local c,y;
  {y,c}=dempiric(x);
 retp(y,CUMSUMC(c));
endp;

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