### Selected R commands and output used in Chapter 3
> data(faithful)
> dim(faithful)
[1] 272   2
> par(mfcol=c(2,2))
> plot(waiting~eruptions,data=faithful,main="linear",pch="*")
> m1 <- lm(waiting~eruptions,data=faithful)
> lines(faithful$eruptions,m1$fitted.values)
> plot(waiting~eruptions,data=faithful,main="bandwidth=0.2",pch="*")
> lines(ksmooth(faithful$eruptions,faithful$waiting,"normal",0.2))
> plot(waiting~eruptions,data=faithful,main="bandwidth=0.5",pch="*")
> lines(ksmooth(faithful$eruptions,faithful$waiting,"normal",0.5))
> plot(waiting~eruptions,data=faithful,main="bandwidth=2",pch="*")
> lines(ksmooth(faithful$eruptions,faithful$waiting,"normal",2))
> par(mfcol=c(1,2))
> require(sm)
> h <- hcv(faithful$eruptions,faithful$waiting,display="lines")
> h
[1] 0.4243375
> sm.regression(faithful$eruptions,faithful$waiting,h=h,xlab="eruptions",
ylab="waiting")

####
> da <- read.table("q-unempdur.txt",header=T)
> y <- da[,2]
> dy <- diff(y)
> ts.plot(dy)
> m1 <- arima(dy,order=c(1,0,0),include.mean=F)
> m1
Call:arima(x = dy, order = c(1, 0, 0), include.mean = F)
Coefficients:
         ar1
      0.5238
s.e.  0.0510

sigma^2 estimated as 0.566:  log likelihood = -314.38,  aic = 632.77
> k1 <- lowess(dy[2:277]~dy[1:276]) # lowess fitting
> require(np)
> Y <- dy[2:277]
> X <- matrix(dy[1:276],276,1)
> m2 <- npreg(txdat=X,tydat=Y,residuals=T)
> names(m2)
 [1] "bw"          "bws"         "pregtype"    "xnames"      "ynames"
 [6] "nobs"        "ndim"        "nord"        "nuno"        "ncon"
[11] "pscaling"    "ptype"       "pckertype"   "pukertype"   "pokertype"
[16] "eval"        "mean"        "merr"        "grad"        "gerr"
[21] "resid"       "ntrain"      "trainiseval" "gradients"   "residuals"
[26] "R2"          "MSE"         "MAE"         "MAPE"        "CORR"
[31] "SIGN"        "rows.omit"   "nobs.omit"   "call"
> m2$bw
[1] 0.4163248
> m2$pckertype
[1] "Second-Order Gaussian"
> fit <- Y-m2$resid
> k2 <- sort(X[,1],index=T)
> tdx <- c(1:278)/4+1948
> par(mfcol=c(2,1))
> plot(tdx,y,xlab='year',ylab='q-dur',type='l')
> plot(tdx[-1],dy,xlab='year',ylab='diff(dur)',type='l')
> par(mfcol=c(1,1))
> plot(X[,1],Y,xlab='Lag-1',ylab='diff(dur)')
> lines(X[,1][k2$ix],fit[k2$ix])
> lines(k1$x,k1$y,lty=2)

####
da <- read.table("SN_y_tot_V2.0.txt",header=T)
sn <- da[,2]; n <- length(sn)
s <- sqrt(var(sn))
y <- sn[3:n]; x <- sn[1:(n-2)]
require(KernSmooth)
par(mfcol=c(2,2))
plot(x,y,xlab="x(t-2)",ylab="x(t)",main=''(a) lowess'')
lines(lowess(x,y),col="red")
m3 <- locpoly(x,y,degree=0,bandwidth=0.5*s)
plot(x,y,xlab="x(t-2)",ylab="x(t)",main=''(b) locpoly-0'')
lines(m3$x,m3$y,col="red")
m4 <- locpoly(x,y,degree=1,bandwidth=0.5*s)
plot(x,y,xlab="x(t-2)",ylab="x(t)",main=''(c) locpoly-1'')
lines(m4$x,m4$y,col=''red'')
m5 <- locpoly(x,y,degree=2,bandwidth=0.5*s)
plot(x,y,xlab="x(t-2)",ylab="x(t)",main=''(d) locpoly-2'')
lines(m5$x,m5$y,col="red")

####
> x <- seq(0,1,0.02); f <- 2*sin(2*pi*x)
> set.seed(11)
> y <- f + rnorm(51)
> require(stats)
> par(mfcol=c(2,2))
> xi <- c(0.32,0.66) # equally-spaced knots
> plot(x,y); abline(v=xi,col="red")
> m1 <- smooth.spline(x,y,df=2)
> lines(m1)
> plot(x,y); abline(v=xi,col="red")
> m2 <- smooth.spline(x,y,df=3)
> lines(m2)
> plot(x,y); abline(v=xi,col="red")
> m3 <- smooth.spline(x,y,df=4)
> lines(m3)
> plot(x,y); abline(v=xi,col="red")
> m4 <- smooth.spline(x,y,df=5)
> lines(m4)

####
> require(stats); require(splines)
> da <- read.table("m-gmsp6708.txt",header=T)
> da1 <- read.table("m-riskfree3m-6708.txt",header=T)
> rf <- da1[,4]/(100*12)
> gm <- da$GM-rf; sp <- da$SP-rf
> plot(sp,gm,xlab="Market excess returns",ylab='GM excess returns')
> m1 <- lm(gm~sp)
> lines(sp,m1$fitted.values)
> m2 <- bs(sp,degree=3) ### B-splines of degree 3.
> m3 <- lm(gm~m2)
> xx <- sort(sp,index=T)
> lines(sp[xx$ix],m3$fitted.values[xx$ix],col=2)

####
da <- read.table("GDPC1.txt",header=T)
gdp <- diff(log(da$rgdp)); n <- length(gdp)
x <- gdp[1:(n-2)]; y <- gdp[3:n]
plot(x,y,xlab='gdp(t-2)',ylab='gdp(t)')
require(splines)
m1 <- smooth.spline(x,y,cv=T)  ## cross-validation
lines(m1)
names(m1)
 [1] "x"        "y"        "w"        "yin"      "data"     "lev"
 [7] "cv.crit"  "pen.crit" "crit"     "df"       "spar"     "lambda"
[13] "iparms"   "fit"      "call"
m1$spar
[1] 1.338832
m2 <- smooth.spline(x,y,df=5)
lines(m2,lty=2) #col="red")
m2$spar
[1] 1.179196
m3 <- smooth.spline(x,y,df=15)
lines(m3,lty=3) #col="blue")
m3$spar
[1] 0.8719036

####
> require(wmtsa)
> plot(ecg)
> m1 <- wavDWT(ecg,n.level=6)
> names(m1)
[1] "data"       "scales"     "series"     "dictionary" "shifted"
[6] "xform"
> plot(m1)
> eda.plot(m1)
> m2 <- wavDWT(ecg,n.level=6,wavelet="haar")
> plot(m2)
> eda.plot(m2)
# Thresholding with "universal" threshold function
> d1 <- wavShrink(m1$data$d1) ## hard threshold
> d1a <- wavShrink(m1$data$d1,shrink.fun="soft")
> par(mfcol=c(3,1))
> plot(m1$data$d1,type="h",ylab="d1 coefficients")
> plot(d1,type="h",ylab='hard')
> plot(d1a,type="h",ylab='soft')
### Reconstruction
> d1 <- wavShrink(m1$data$d1)  ## Hard thresholding
> d2 <- wavShrink(m1$data$d2)
> d3 <- wavShrink(m1$data$d3)
> d4 <- wavShrink(m1$data$d4)
> m1copy <- m1
> m1copy$data$d1 <- d1
> m1copy$data$d2 <- d2
> m1copy$data$d3 <- d3
> m1copy$data$d4 <- d4
> ecg1 <- reconstruct(m1copy)
> ts.plot(ecg1,ylab='ecg',col="red")
> points(1:2048,as.numeric(ecg),cex=0.4)

####
> require(mda)
> m3.bruto <- bruto(X[,-1],X[,1])
> m3.bruto$type
 pcgasm1  pcgasm2    pcoil  pcoilm1
  smooth excluded   smooth   linear
Levels: excluded linear smooth

####
require(gam)
da <- read.csv(``w-gasoil-pc.csv'')
X <- data.frame(da[,-1])
m2.lm <- lm(pcgas~-1+pcgasm1+pcgasm2+pcoil+pcoilm1,data=X)
m2 <- gam(pcgas~-1+s(pcgasm1,4)+s(pcgasm2,4)+s(pcoil,4)+s(pcoilm1,4),data=X)
summary(m2)
m2.gam <- gam(pcgas~-1+s(pcgasm1,4)+s(pcoil,4)+s(pcoilm1,4),data=X)
par(mfcol=c(2,2))
plot(m2.gas)
par(mfcol=c(1,1))
anova(m2.lm,m2.gam)
require(mda)
m3.bruto <- bruto(X[,-1],X[,1])
m3.bruto$type

####
> require(dr)
> da <- read.csv("w-gasoil-pc.csv")
> head(da)
  N  pcgas pcgasm1 pcgasm2 pcoil pcoilm1
1 1 -0.038  -0.050  -0.047 -2.10   -3.25
....
> X <- da[,-1]
> m1 <- dr(pcgas~.,data=data.frame(X))
> summary(m1)
Call: dr(formula = pcgas ~ ., data = data.frame(X))

Method:  sir with 8 slices, n = 478.
Slice Sizes:
59 59 59 60 61 59 61 60

Estimated Basis Vectors for Central Subspace:
            Dir1      Dir2     Dir3      Dir4
pcgasm1 0.996244 -0.037141  0.65816  0.902187
pcgasm2 0.083244  0.999245 -0.75199 -0.431118
pcoil   0.023252 -0.011099 -0.02416 -0.002282
pcoilm1 0.005176 -0.002633  0.02744 -0.013846

              Dir1    Dir2     Dir3     Dir4
Eigenvalues 0.4983 0.03988 0.004562 0.001055
R^2(OLS|dr) 0.9975 0.99855 0.999395 1.000000

Large-sample Marginal Dimension Tests:
                Stat df p.value
0D vs >= 1D 259.9249 28  0.0000
1D vs >= 2D  21.7489 18  0.2433
2D vs >= 3D   2.6851 10  0.9879
3D vs >= 4D   0.5043  4  0.9731

####
### Further Analysis
> TZ <- as.matrix(X[,-1])%*%m1$evectors
> idx <- c(1:478)[TZ[,1] < -0.2]
> R1 <- rep(0,478)
> R1[idx] <- 1
> LTZ <- TZ[,1]*R1
> m2 <- lm(X$pcgas~-1+TZ[,1]+LTZ)
> summary(m2)
Call: lm(formula = X$pcgas ~ -1 + TZ[, 1] + LTZ)

Coefficients:
        Estimate Std. Error t value Pr(>|t|)
TZ[, 1]  0.36633    0.02017  18.165  < 2e-16 ***
LTZ      0.09613    0.03345   2.874  0.00424 **
---
Residual standard error: 0.04152 on 476 degrees of freedom
Multiple R-squared:  0.5697,    Adjusted R-squared:  0.5679

> Box.test(m2$residuals,lag=12,type='Ljung')
        Box-Ljung test
data:  m2$residuals
X-squared = 10.353, df = 12, p-value = 0.585

####