### R code from vignette source 'C:/Users/Marie/Documents/DSP/cviceni/M5201/M5201cv04.Rnw'
### Encoding: CP1250

###################################################
### code chunk number 1: M5201cv04.Rnw:55-56
###################################################
options(continue=" ")


###################################################
### code chunk number 2: M5201cv04.Rnw:59-61
###################################################
lattice.options(default.theme = "col.whitebg")
options(SweaveHooks = list(fig = function() par(mfrow=c(1,1))),width=80)


###################################################
### code chunk number 3: M5201cv04.Rnw:399-401
###################################################
fileDat<-paste(data.library,"treerings.txt",sep="")
(TXT<-paste(scan(fileDat,what="",nlines=1),collapse=" "))


###################################################
### code chunk number 4: M5201cv04.Rnw:408-410
###################################################
rings<-read.table(fileDat,header=F,skip=2)
str(rings)


###################################################
### code chunk number 5: M5201cv04.Rnw:417-420
###################################################
getOption("SweaveHooks")[["fig"]]()
ringsTS<-ts(rings,start=1753,frequency=1)
par(mar=c(2,2,1,0)+0.5)
plot(ringsTS,main=TXT,cex.main=0.85)


###################################################
### code chunk number 6: M5201cv04.Rnw:434-441
###################################################
y <- as.vector(ringsTS)
x <- as.vector(time(ringsTS))
n<-length(y)
nn<-300
data<-data.frame(x,y)
LinTrend <- lm( y ~ x,data=data)
print(summary(LinTrend))


###################################################
### code chunk number 7: M5201cv04.Rnw:450-458
###################################################
getOption("SweaveHooks")[["fig"]]()
new <- data.frame(x = seq(x[1],x[n],length.out=nn))
pred.w.plim <- predict(LinTrend, new, interval="prediction")
pred.w.clim <- predict(LinTrend, new, interval="confidence")
par(mar=c(2,2,1,0)+0.5)
plot(x,y,type="l")
matlines(new$x,cbind(pred.w.clim, pred.w.plim[,-1]),col=c(2,3,3,4,4),
        lty=c(1,2,2,3,3), type="l", ylab="predicted y")
title(main=TXT,cex.main=0.85)


###################################################
### code chunk number 8: M5201cv04.Rnw:473-474
###################################################
res<-resid(LinTrend)


###################################################
### code chunk number 9: M5201cv04.Rnw:482-484
###################################################
fileFun<-paste(data.library,"FunkceM5201.R",sep="")
source(fileFun)


###################################################
### code chunk number 10: M5201cv04.Rnw:492-494
###################################################
getOption("SweaveHooks")[["fig"]]()
pdgV<-tperiod(res,fig=TRUE)
str(pdgV)


###################################################
### code chunk number 11: M5201cv04.Rnw:515-517
###################################################
(TkF<-n/pdgV$kF)
(TkS<-n/pdgV$kS)


###################################################
### code chunk number 12: M5201cv04.Rnw:543-545
###################################################
TrigoTrend<-trigofit(res)
str(TrigoTrend)


###################################################
### code chunk number 13: M5201cv04.Rnw:553-554
###################################################
FisherPeriods<-2*pi/TrigoTrend$signFreq


###################################################
### code chunk number 14: M5201cv04.Rnw:574-580
###################################################
getOption("SweaveHooks")[["fig"]]()

TTIME<-seq(1,length(res),length=nn)
predTrigoTrend<-trigoval(TTIME,A=TrigoTrend$signA,B=TrigoTrend$signB,freq=TrigoTrend$signFreq)
par(mar=c(2,2,1,0)+0.5)
plot(res~x, type="l",main=TXT,cex.main=0.85)
lines(new$x,predTrigoTrend,col="dodgerblue",lwd=2)


###################################################
### code chunk number 15: M5201cv04.Rnw:595-597
###################################################
betaD<-damsleth(res,pstlevel=0.05)
str(betaD)


###################################################
### code chunk number 16: M5201cv04.Rnw:608-609
###################################################
(DamslethPeriods<-2*pi/betaD$freq)


###################################################
### code chunk number 17: M5201cv04.Rnw:618-629
###################################################
getOption("SweaveHooks")[["fig"]]()
predTrigD<-trigoval(TTIME,A=betaD$A,B=betaD$B,freq=betaD$freq)
par(mar=c(2,2,1,0)+0.5)
plot(res~x, type="l",main=TXT,cex.main=0.85)
lines(new$x,predTrigoTrend,col="dodgerblue",lwd=2)
lines(new$x,predTrigD,lwd=2,col="red",lty=2)
TFP<-paste("Fisher's periods:      ",
     paste(round(FisherPeriods,1),collapse=",  "))
TDS<-paste("Damsleth-Spjotvoll:",
     paste(round(DamslethPeriods,1),collapse=",  "))
legend("topright", legend = c(TFP,TDS),bty="n",
                    lty = c(1,2), xjust = 1, yjust = 1,col=c("dodgerblue","red"))


###################################################
### code chunk number 18: M5201cv04.Rnw:643-651
###################################################
getOption("SweaveHooks")[["fig"]]()
pred.lintrend<-predict(LinTrend, new)
par(mar=c(2,2,1,0)+0.5)
plot(y~x, type="l",main=TXT,cex.main=0.85)
lines(new$x,predTrigoTrend+pred.lintrend,col="dodgerblue",lwd=2)
lines(new$x,predTrigD+pred.lintrend,lwd=2,col="red",lty=2)
lines(new$x,pred.lintrend,lwd=1.5,col="darkgreen",lty=3)
legend("topright", legend = c(TFP,TDS,"linear trend"),bty="n",
                    lty = c(1,2,3), xjust = 1, yjust = 1,col=c("dodgerblue","red","darkgreen"))


