calendar 1881 1 12
allocate 2003:4
open data pe-sp500.xls
data(format=xls, org=col) / pe10

com end = 2003:4
com start = 1881:1
set time = t
sta(noprint) time  ; com iloops = 200

set y start end = pe10
set pe start end = pe10

dif y / dy

com rss_min = 10000000000.
do j = 1,iloops
       com k = j/20.
       sta(noprint) y ; com nobs = %nobs-1
       set sin1 = sin((2*%pi/nobs)*k*time)
       set cos1 = cos((2*%pi/nobs)*k*time)

		com lags = 13, hit = 0
 		until hit.eq.1 {
		com lags = lags-1
		lin(noprint) dy; # dy{1 to lags} sin1 cos1 constant y{1}
		if abs(%tstats(lags)).ge.(1.95) ; com hit = 1
		}
        if lags>0 {
	   lin(noprint) dy ; # dy{1 to lags} sin1 cos1 constant y{1}  ; }
        if lags.le.0 {
           lin(noprint) dy ; # sin1 cos1 constant y{1}  ; }
	if %rss < rss_min
	    {
	    compute kmin = k
            compute rss_min = %rss
            }
end do j
dis kmin
  set sin1  = sin((2*%pi/nobs)*kmin*time)
  set cos1  = cos((2*%pi/nobs)*kmin*time)

  	com lags = 13, hit = 0
	until hit.eq.1 {
	com lags = lags-1
	lin(noprint) dy ; # dy{1 to lags} sin1 cos1 constant y{1}
	if abs(%tstats(lags)).ge.(1.95) ; com hit = 1
	}
         end until

*  if lags>0 {
  	linreg dy ;# constant sin1 cos1 y{1} dy{1 to lags}
*  if lags<=0 {
*  linreg(noprint) dy ; # constant sin1 cos1 y{1}
  exc(noprint); # sin1 cos1
  com tau = %tstats(4)


linreg(noprint) dy ; # constant y{1} dy{1 to lags}
com df = %tstats(2)
dis  '%cdstat tau kmin lags df'
dis  %cdstat tau kmin lags df



  linreg y ;# constant sin1 cos1 y{1} y{1 to lags+1}
  summarize(noprint) ; # y{1 to lags+1} ; com dem = 1-%sumlc
  set Fourier start end = (%beta(1) + %beta(2)*sin1 + %beta(3)*cos1)/dem

*end do start

label pe
# 'PE'

*label r3 tbill r1
*# '3-year' 'T-bill' '1-year'


*gra(Hea="Figure 1: .. ",key=below,nokbo,patt,vla='percent per year') 3 ; # r3 1991:1 end; # r1 1991:1 end; # tbill 1991:1 end
* grparm header 22

*set spread = r3 - r1

gra(hea="PE and the Fourier Term", $
	patterns,vla="PE",key=below,nokbo) $
	2 ; # PE 1881:1 *  ; # Fourier 1881:1 *

gra(hea=" ", $
	patterns,vla="PE",key=below,nokbo) $
	2 ; # PE 1881:1 *  ; # Fourier 1881:1 *


/*
set y = pe10; dif y / dy
set sin1  = sin((2*%pi/nobs)*time)
set cos1  = cos((2*%pi/nobs)*time)

lin dy ; # y{1} constant sin1 cos1 dy{1 to 11}
lin(noprint) y ; # constant sin1 cos1 y{1 to 12}
summarize(noprint) # y{1 to 12}
set y = pe


*lin y ; # constant y{1 to 2}
*source c:\winrats\baiperron.src
*@baiperron(maxbreaks=5,minspan=12) y
*# constant y{1 to 2}

*dis %datelabel(59) %datelabel(132) %datelabel(150)
*set c1 = %if(t<59,1.,0.)
*set c2 = %if(t.ge.59.and.t<132,1.,0.)
*set c3 = %if(t.ge.132.and.t<150,1.,0.)
*set c4 = %if(t.ge.150,1.,0.)
*lin y ; # c1 c2 c3 c4 y{1 to 2}
*set bp_int = (%beta(1)*c1+%beta(2)*c2+%beta(3)*c3+%beta(4)*c4)/(1-%beta(5)-%beta(6))

*label bp_int intercept
*# 'Bai-Perron' 'Fourier'

*gra(hea="Figure 5: The R3 - R1 Spread and the Two Intercepts", $
*	patterns,vla="Spread",key=below,nokbo,noaxis) $
*	3 ; # spread 1991:1 *  ; # intercept 1991:1 * 1 ; # bp_int 1991:1 *

*/