/*
This is a modified GAUSS program file based on the code originally written by

Bruce E. Hansen
Department of Economics
Social Science Building
University of Wisconsin
Madison, WI 53706-1393
bhansen@ssc.wisc.edu
http://www.ssc.wisc.edu/~bhansen/

replicating the estimation, testing and graphs reported in
"Threshold Effects in Non-Dynamic Panels:
Estimation, Testing and Inference" 

-----------------------------------------------------------------

This GAUSS program file replicates the estimation, testing and graphs reported in
"Inflation Thresholds and Relative Price Variability: 
Evidence from U.S. Cities"

The code is augmented by regime intercepts in all but one regime. Hansen (1999) 
allows for more than one regressor. This program was modified for the case of
a single regressor as in our application. 

For questions regarding these extensions, please contact

Alexander Bick
Department of Economics
Postbox 28
Goethe-University Frankfurt
House of Finance
Grüneburgplatz 1
60323 Frankfurt am Main
bick@wiwi.uni-frankfurt.de

The bootstrap procedure has also been modified. Hansen (1999) groups the 
regression residuals by individual e_i={e_i1, e_i2,..., e_iT} and takes the 
sample {e_1, e_2,..., e_N} with size N as the empirical distribution. Since N is limited 
but T is large in our empirical analysis, we treat the 
sample {e_11, e_12,..., e_1T,..., e_i1,..., e_iT,..., e_N1,...,e_NT} as empirical distribution 
to be used for bootstrapping. For the bootstrap procedure, the variable x_it
and the threshold variable q_it are given, i.e. their values are fixed in repeated bootstrap samples. We
take with replacement a sample of size NT from the empirical distribution and create a bootstrap
sample under the null hypothesis of no threshold. This bootstrap sample is used to estimate the
model under H0 and H1 and to calculate the bootstrap value of the likelihood ratio statistic F1. 
This procedure is frequently repeated - 1000 bootstrap replications in our application - and
the bootstrap estimate of the asymptotic p-value for F1 under H0 is the percentage of draws for
which the simulated likelihod ratio statistic exceeds the actual statistic. The null hypothesis of no
threshold effect is rejected if the p-value is smaller than the desired significance level.
This procedure has been implemented by Juliane Scharff. 
*/

load rpv_inf[560,2]=rpv_inf.txt;
inf         =rpv_inf[.,1];               @ inflation                             @
rpv         =rpv_inf[.,2];               @ rpv                                   @

reg_tresh   =inf;               	 @ regressor for tresholdvariable        @
tresh       =inf;                  	 @ tresholdvariable                      @


t = 40;             
nt= rows(rpv);      
n = nt/t;           


conf_lev = .95;      @ confidence level for threshold                @
_vgraph = 1;         @ set to 1 to graph likelihood ratios           @
_boot_1 = 1000;      @ # of replications, 0 for no bootstrap, single (300) @
_boot_2 = 1000;      @ # of replications, 0 for no bootstrap, double (300) @
_boot_3 = 1000;      @ # of replications, 0 for no bootstrap, triple (300) @
_trim_1 = .05;       @ percentage to trim before search, single      @
_trim_2 = .05;       @ percentage to trim before search, double      @
_trim_3 = .05;       @ percentage to trim before search, triple      @

output file=rpv_inf.out reset;  output off;

max_lag = 0;        
tt = t-max_lag;     
ty = n*(t-max_lag-1);

y  = rpv;  
yt = tr(y);
cf = reg_tresh;
ct = tr(cf);

d1 = tresh;      @ set to threshold variable @

thresh = d1;            
dd = unique(thresh,1);  
qn = rows(dd);         
qnt1 = ceil(qn*_trim_1);      
                                
qn_help=qn-qnt1;
qq1 = dd[qnt1:qn_help];   
qn1 = rows(qq1);      
cc = -2*ln(1-sqrt(conf_lev));   

output on;
"Output:                headline_only";

"Number of Cities        " n;
"Number of Months used   " tt;
"Total Observations      " ty;
"Number of Quantiles     " qn;
"Confidence Level        " conf_lev;
"";"";
"*******************************************************";
"";"";
sse0 = sse_calc(yt,ct);
"Zero Threshold Model";
"Sum of Squared Errors                   " sse0;
"";"";
"*******************************************************";
"";"";

"Single Threshold Model";
"";
output off;
rhat1 = model(0,_trim_1,_boot_1,0);
output on;
"*******************************************************";
"";"";

"Double Threshold Model";
"Trimming Percentage    " _trim_2;
"";
"First Iteration";
output off;
rhat2 = model(rhat1,_trim_2,_boot_2,2);
output on;
"Second Iteration";
output off;
rhat1 = model(rhat2,_trim_2,0,1);
output on;
"";"";
"*******************************************************";
"";"";

"Triple Threshold Model";
"Trimming Percentage    " _trim_3;
"";
output off;
rhat3 = model(rhat1|rhat2,_trim_3,_boot_3,3);
output on;
"";"";
"*******************************************************";
"";"";
output off;


/****************** PROCS ************************/

proc tr(y);
 local yfm,yf;
 yf = reshape(y,n,tt);
 yfm = yf - meanc(yf');
 yfm = yfm[.,1:tt-1];
retp(vec(yfm'));
endp;

proc lag_v(x,lagn);
local y;
  y = reshape(x,n,t);
  y = y[.,1+max_lag-lagn:t-lagn];
retp(vec(y'));
endp;

proc sse_calc(y,x);
local e;
  e = y - x*(y/x);
retp(e'e);
endp;


proc thr_sse(y,q,r);
local n,sse,qi,rr,d,xx,dum_s,j;
n = rows(q);
sse = zeros(n,1);
qi = 1; do while qi<=n;
  if r==0; rr = q[qi]; else; rr = r|q[qi]; endif;
  rr = sortc(rr,1);
  xx = ct;
  j = 1; do while j <= rows(rr);
     d = (thresh .< rr[j]);
     xx = xx~tr(cf.*d)~tr(d);
  j=j+1;endo;
  sse[qi] = sse_calc(y,xx);
qi = qi+1;endo;
retp(sse);
endp;


proc (2) = r_est(y,r,_trim);
local qq,rr,i,nn,ii,sse,rihat,qnt;
if maxc(r) .== 0;
  qq = qq1; rr = 0;
else;
  rr = sortc(r,1);
  i = seqa(1,1,qn1);
  nn = sumc(qq1 .< (rr'))';
  qnt = qn*_trim;
  ii = ((i .<= (nn+qnt))-(i .<= (nn-qnt)))*ones(rows(rr),1);
  qq = delif(qq1,ii);
endif;
sse = thr_sse(y,qq,rr);
rihat = minindc(sse);
retp(sse[rihat],qq[rihat]);
endp;


proc (1) = model(r,_trim,rep,it);
local qq,rr,i,nn,ii,sse,rihat,rhat,sse1,lr,rhats,xx,crits,
rrr,e,sse0,yp,lrt,stats,j,eb,yb,yp_b,lrt_b,jj,dd,d,xxi,beta,
sehet,sehomo,titname,tit,xname,xlab,nr,rhat_b,qnt,
z; 

if maxc(r)==0;
  qq = qq1; rr = 0;
else;
  rr = sortc(r,1);
  i = seqa(1,1,qn1);
  nn = sumc(qq1 .< (rr'));
  qnt = qn*_trim;
  ii=((i.<=(nn+qnt)')-(i.<=(nn-qnt)'))*ones(rows(rr),1);
  qq = delif(qq1,ii);
endif;
sse = thr_sse(yt,qq,rr);
rihat = minindc(sse);
rhat = qq[rihat];
sse1 = sse[rihat];
lr = (sse/sse1 - 1)*ty;
rhats = selif(qq,(lr .< cc));
if _vgraph==1;
  library pgraph; 
  graphset;
  if     it==0; 
    titname="Figure 1\LConfidence Interval Construction in Single Threshold Model";
    xname="Threshold Parameter";
  elseif it==1; 
    titname = "Figure 3\LConfidence Interval Construction in Double Threshold Model"; 
    xname="First Threshold Parameter";
  elseif it==2; 
    titname = "Figure 2\LConfidence Interval Construction in Double Threshold Model"; 
    xname="Second Threshold Parameter";
  elseif it==3; 
    titname = "Confidence Interval Construction in Triple Threshold Model"; 
    xname="Third Threshold Parameter";
  endif;
  
  ylabel("Likelihood Ratio");
  title("Testtitel");
  xlabel("x-Variable");
  xy(qq,lr~(ones(rows(qq),1)*cc)); 
endif;
output on;
if maxc(r) .ne 0;
  "Fixed Thresholds       " rr';
  rrr = sortc((rr|rhat),1);
else;
  rrr = rhat;
endif;
"Threshold Estimate     " rhat;
"Confidence Region      " minc(rhats)~maxc(rhats);
"Trimming Percentage    " _trim;
"";"";

nr = rows(rrr);
xx = ones(rows(yt),1);
dd = zeros(rows(thresh),nr);
j=1; do while j<=nr;
  dd[.,j] = (thresh .< rrr[j]);
  d = dd[.,j];
    if j>1;
    d = d - dd[.,j-1];
  endif;
    xx = xx~tr(cf.*d)~tr(d); @xx = tr(cf.*d);@ 
j=j+1;endo;
d = 1-dd[.,nr];
xx = xx~tr(cf.*d);


z = cols(xx);
xx  = xx[.,2:z];
xxi = invpd(moment(xx,0));
beta = xxi*(xx'yt);
e = yt - xx*beta;
sehet = sqrt(diag(xxi*moment(xx.*e,0)*xxi));
sehomo = sqrt(diag(xxi*(e'e)/(ty-n-cols(xx))));
beta = beta~sehomo~sehet;
"Thresholds";
rrr';
"";
"Regime-independent Coefficients, standard errors, het standard errors";
@beta[1:k,.];@
"";
"Regime-dependent Coefficients, standard errors, het standard errors";
beta;
@beta[k+1:k+nr+1,.];@
"";"";
output off;
if rep .> 0;
  xx = ct;
  if maxc(rr) .ne 0;
    j = 1; do while j <= rows(rr);
     d = (thresh .< rr[j]);
     xx = xx~tr(cf.*d)~tr(d);
    j=j+1;endo;
  endif;
  yp = xx*(yt/xx);
  e = yt-yp;
  sse0 = e'e;
  lrt = (sse0/sse1-1)*ty;
  output on;
  "LR Test for threshold effect  " lrt;
  output off;
  "";"";
  stats = zeros(rep,1);
  j=1; do while j<=rep;
    eb   = e;
    yb   = yp + eb[ceil(rndu(ty,1)*ty)];
    sse0 = sse_calc(yb,ct);
    {sse1,rhat_b} = r_est(yb,0,_trim);
    rrr = rhat_b;
    if maxc(r) .ne 0;
      jj = 1; do while jj<=rows(r);
        sse0 = sse1;
        {sse1,rhat_b} = r_est(yb,rrr,_trim);
        rrr = rrr|rhat_b;
      jj = jj + 1; endo;
    endif;
    lrt_b = (sse0/sse1-1)*ty;
    stats[j] = lrt_b;
    "Bootstrap Replication " j~lrt_b;
  j=j+1;endo;
  "";"";
  stats = sortc(stats,1);
  crits = stats[ceil((.90|.95|.99)*rep)];
  output on;
  "Number of Bootstrap replications   " rep;
  "Bootstrap p-value                  " meanc(stats .> lrt);
  "Critical Values   " crits;
  "";"";
  output off;
endif;
retp(rhat);
endp;


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