# But let's try a different tack.  The heart and respiration cycles are roughly\
# sinusoidal, but the weird cycle is not.  So let's revert to modelling just\
# series #1 but use more harmonics in the weird cycle.  That way the evolution\
# variance need not be so large and maybe won't pick up the signal.\
\
build.cyc5 <- function ( par ) \{\
  # build a model with four components.  One is locally linear.  The\
  # other three are cyclic --- for respiration, heartbeat, and weird.\
  # par[1] is the log of the observation variance\
  # par[2] is the log of the locally linear innovation variance\
  # par[3] is the log of the heartbeat innovation variance\
  # par[4] is the log of the heartbeat innovation variance\
  # par[5] is the log of the weird innovation variance\
  return (   dlmModPoly ( order=2, dV=exp(par[1]), dW=c(0,exp(par[2])) )\
           + dlmModTrig ( tau=2.78,  q=1, dV=0, dW=exp(par[3]) )\
           + dlmModTrig ( tau=18.75, q=1, dV=0, dW=exp(par[4]) )\
           + dlmModTrig ( tau=117,   q=2, dV=0, dW=exp(par[5]) )\
         )\
\}\
system.time ( mle.cyc5 <- dlmMLE ( A535$Region.1, c(0,0,0,0,0), build.cyc5 ) )\
exp ( mle.cyc5$par )\
exp ( mle.cyc3$par )\
exp ( mle.cyc2$par )\
exp ( mle.cyc1$par )\
\
mod.cyc5 <- build.cyc5 ( mle.cyc5$par )\
filt.cyc5 <- dlmFilter ( A535$Region.1, mod.cyc5 )\
smoo.cyc5 <- dlmSmooth ( filt.cyc5 )\
\
par ( mfrow=c(3,2) )\
plot.smoo.component ( smoo.cyc5, c(1,0,1,0,1,0,1,0,1,0), A535$Region.1,\
                      plot.data=T, pch=1, xlab="time", ylab="", \
                      main="Fit")\
plot.smoo.component ( smoo.cyc5, c(1,0,0,0,0,0,0,0,0,0), A535$Region.1, \
                      plot.data=T, pch=1, xlab="time", ylab="", \
                      main="Level")                      \
plot.smoo.component ( smoo.cyc5, c(0,1,0,0,0,0,0,0,0,0), A535$Region.1, \
                      pch=1, xlab="time", ylab="", main="Slope")\
abline ( h=0, lty=2 )\
plot.smoo.component ( smoo.cyc5, c(0,0,1,0,0,0,0,0,0,0), A535$Region.1, \
                      pch=1, xlab="time", ylab="", main="Heartbeat")\
plot.smoo.component ( smoo.cyc5, c(0,0,0,0,1,0,0,0,0,0), A535$Region.1, \
                      pch=1, xlab="time", ylab="", main="Respiration")\
plot.smoo.component ( smoo.cyc5, c(0,0,0,0,0,0,1,0,1,0), A535$Region.1, \
                      pch=1, xlab="time", ylab="", main="Weird")\
# That looks much better.  Check residuals\
resids.raw <- residuals ( filt.cyc5, type="raw", sd=F )\
resids.std <- residuals ( filt.cyc5, type="st", sd=F )\
\
par ( mfrow=c(3,2) )\
plot ( resids.raw, type="l" )\
plot ( resids.std, type="l" )\
acf ( resids.std )\
acf ( resids.std, lag.max=200 )\
qqnorm ( resids.std )\
\
# There is still some cyclic behavior with a period of about 10.  That's very\
# close to half the respiration period.  I'm going to add a harmonic to the\
# respiration\
build.cyc6 <- function ( par ) \{\
  # build a model with four components.  One is locally linear.  The\
  # other three are cyclic --- for respiration, heartbeat, and weird.\
  # par[1] is the log of the observation variance\
  # par[2] is the log of the locally linear innovation variance\
  # par[3] is the log of the heartbeat innovation variance\
  # par[4] is the log of the heartbeat innovation variance\
  # par[5] is the log of the weird innovation variance\
  return (   dlmModPoly ( order=2, dV=exp(par[1]), dW=c(0,exp(par[2])) )\
           + dlmModTrig ( tau=2.78,  q=1, dV=0, dW=exp(par[3]) )\
           + dlmModTrig ( tau=18.75, q=2, dV=0, dW=exp(par[4]) )\
           + dlmModTrig ( tau=117,   q=2, dV=0, dW=exp(par[5]) )\
         )\
\}\
mle.cyc6 <- dlmMLE ( A535$Region.1, c(0,0,0,0,0), build.cyc6 )\
exp ( mle.cyc6$par )\
exp ( mle.cyc5$par )\
\
mod.cyc6 <- build.cyc6 ( mle.cyc6$par )\
filt.cyc6 <- dlmFilter ( A535$Region.1, mod.cyc6 )\
smoo.cyc6 <- dlmSmooth ( filt.cyc6 )\
\
par ( mfrow=c(3,2) )\
plot.smoo.component ( smoo.cyc6, c(1,0,1,0,1,0,1,0,1,0,1,0), A535$Region.1,\
                      plot.data=T, pch=1, xlab="time", ylab="", \
                      main="Fit")\
plot.smoo.component ( smoo.cyc6, c(1,0,0,0,0,0,0,0,0,0,0,0), A535$Region.1, \
                      plot.data=T, pch=1, xlab="time", ylab="", \
                      main="Level")                      \
plot.smoo.component ( smoo.cyc6, c(0,1,0,0,0,0,0,0,0,0,0,0), A535$Region.1, \
                      pch=1, xlab="time", ylab="", main="Slope")\
abline ( h=0, lty=2 )\
plot.smoo.component ( smoo.cyc6, c(0,0,1,0,0,0,0,0,0,0,0,0), A535$Region.1, \
                      pch=1, xlab="time", ylab="", main="Heartbeat")\
plot.smoo.component ( smoo.cyc6, c(0,0,0,0,1,0,1,0,0,0,0,0), A535$Region.1, \
                      pch=1, xlab="time", ylab="", main="Respiration")\
plot.smoo.component ( smoo.cyc6, c(0,0,0,0,0,0,0,0,1,0,1,0), A535$Region.1, \
                      pch=1, xlab="time", ylab="", main="Weird")\
# Best so far.  Check residuals\
resids.raw <- residuals ( filt.cyc6, type="raw", sd=F )\
resids.std <- residuals ( filt.cyc6, type="st", sd=F )\
\
par ( mfrow=c(2,2) )\
plot ( resids.raw, type="l" )\
plot ( resids.std, type="l" )\
acf ( resids.std, lag.max=200 )\
qqnorm ( resids.std )\
# Not perfect.  But I like this model\
}