Draw the graph - Actual, Theory and Residual.

plot(merge(merge(apply.quarterly(SP5[k2k],mean)[,1],predict(lm(apply.quarterly(SP5[k2k],mean)[,1] ~ PAq[k2k] * UCq[k2k] * G[k2k]*CSq[k2k] - UCq[k2k] -G[k2k] - PAq[k2k]*G[k2k] - UCq[k2k]*G[k2k]*CSq[k2k]))),
+ residuals(lm(apply.quarterly(SP5[k2k],mean)[,1] ~ PAq[k2k] * UCq[k2k] * G[k2k]*CSq[k2k] - UCq[k2k] -G[k2k] - PAq[k2k]*G[k2k] - UCq[k2k]*G[k2k]*CSq[k2k]))))
axis(2,at=c(100,50,0,-50,-100,-150))
axis(4,at=c(100,50,0,-50,-100,-150))

black = actual
red = theory
green = residual

# or draw horizontal line at the designated x-value, please not that
# the index must be synchronized with the main graph, which is equal to
#  seq(as.Date("2000-03-01"),as.Date("2017-12-01"),by='quarters'))
# no to k2k.

plot(merge(merge(apply.quarterly(SP5[k2k],mean)[,1],predict(lm(apply.quarterly(SP5[k2k],mean)[,1] ~ PAq[k2k] * UCq[k2k] * G[k2k]*CSq[k2k] - UCq[k2k] -G[k2k] - PAq[k2k]*G[k2k] - UCq[k2k]*G[k2k]*CSq[k2k]))),
+            + residuals(lm(apply.quarterly(SP5[k2k],mean)[,1] ~ PAq[k2k] * UCq[k2k] * G[k2k]*CSq[k2k] - UCq[k2k] -G[k2k] - PAq[k2k]*G[k2k] - UCq[k2k]*G[k2k]*CSq[k2k]))),ylim=c(-200,3000))
# addSeries(as.xts(rep(-100,72),index(apply.quarterly(SP5[k2k],mean)[,1])),on=1,ylim=c(-200,2700),lwd=1)
addSeries(as.xts(rep(-100,length(index(apply.quarterly(SP5[k2k],mean)[,1]))),index(apply.quarterly(SP5[k2k],mean)[,1])),on=1,ylim=c(-200,3000),lwd=1)

S&P 500 for next 9 qtrs with and w/o case shiller price index

kg <- "1992-01-01::2017-12-31"
kikan <- "1992-01-01::2017-12-31"
k2k <- "2000-01-01::2017-12-31"   #GPUC model is based data only from 2000 because of cs index availability
# prepare parameters.
l <- 48  # # of months to predict
r <- 1.05 # pesumed GDP growth rate
# get data from FRED
getSymbols("GDP",src="FRED")
getSymbols("PAYEMS",src="FRED")
getSymbols("UNDCONTSA",src="FRED")
getSymbols('SPCS10RSA',src='FRED')
CS <- SPCS10RSA
PA <- PAYEMS
UC <- UNDCONTSA
# download data from yahoo finace before this.
SP5 <- as.xts(read.zoo(read.csv("~/SP5.csv")))

G <- GDP
i <- seq(2,l/3,1)  # seq of quarters to predict
d <- as.Date(as.yearqtr(seq(Sys.Date(),as.Date("2100-12-31"),by="quarters")[i])) # pick up the first day of each quarters.
# for the case the current time is more than 6 months away from the last GDP data, the below is to adjust # of months for GDP forecast.
### caution still in debug
if(floor(MonthsBetween(mondate(last(index(G))),mondate(Sys.Date()))) > 5){ i <- i+1}
### caution ends
m_GDP <- as.xts(as.vector(last(GDP)) * r**(i/4),d)
m_PA <- (as.xts(forecast(auto.arima(PA),h=l)$mean[1:l],as.Date(as.yearmon(seq(mondate(index(last(PA)))+1,by=1,length.out=l))))[(3-month(index(last(PA))) %% 3) + seq(1,l-3,3)])[d]
m_UC <- (as.xts(forecast(auto.arima(UC),h=l)$mean[1:l],as.Date(as.yearmon(seq(mondate(index(last(UC)))+1,by=1,length.out=l))))[(3-month(index(last(UC))) %% 3) + seq(1,l-3,3)])[d]
m_CS <-
(as.xts(forecast(auto.arima(CS),h=l)$mean[1:l],as.Date(as.yearmon(seq(mondate(index(last(CS)))+1,by=1,length.out=l))))[(3-month(index(last(CS))) %% 3) + seq(1,l-3,3)])[d]

# predict next 9 qtrs. by GPU model
my_sp5(kikan,m_GDP[1:9],m_PA_backup[1:9],m_UC[1:9])
# by GPUC model borrow index from PA for case shiller index.
# use "as.xts()", not "xts()"
# this assumes case shiller index of the next qtr is 223 and it will add 2 pts. each qtr.
my_sp5cs(k2k,m_GDP[1:9],m_PA[1:9],m_UC[1:9],as.xts(seq(223,223+2*8,2),index(m_PA[1:9])))
# or use m_CS[1:9]

draw the graph to record blood pressure vol.2

Please refer to 2018/2/20 for the detail


#
# read data as csv format and convert to xts to plot
#
Sys.setenv(TZ=Sys.timezone())
bp <- read.csv("~/Downloads/bp - シート1.csv")
system("rm \"$HOME/Downloads/bp - シート1.csv\"")
bp.xts <- xts(bp[,c(-1,-2)],as.POSIXct(paste(bp$Date,bp$Time,sep=" "),tz=Sys.timezone()),tz=Sys.timezone())
# bp.xts <- xts(bp[,c(-1,-2)],as.POSIXct(paste(bp$Date,bp$Time,sep=" ")))
# bp.xts <- xts(bp[,c(-1,-2)],as.POSIXct(paste(bp$Date,bp$Time,sep=" ")),tz="UTC")
# weekly average
apply.weekly(bp.xts[bp.xts$High > 95],mean)
# make UTC based weekly data
# apply.weekly(as.xts(bp[,c(-1,-2)],as.POSIXct(paste(bp$Date,bp$Time,sep=" "),tz="UTC")),mean)
#  draw the graph
plot(bp.xts[,-3][bp.xts$High > 95],col = c("red", "blue"),lwd=c(3,3,2,2),major.ticks='days',grid.ticks.on='days',type='p',ylim=c(60,160))
# draw a horizontal line at 130 as the benchmark
addSeries(xts(rep(130,length(index(bp.xts[bp.xts$High > 95]))),index(bp.xts[bp.xts$High > 95])),ylim=c(60,160),on=1,col=5,lwd=1)
# draw another line at 85.
addSeries(xts(rep(85,length(index(bp.xts[bp.xts$High > 95]))),index(bp.xts[bp.xts$High > 95])),ylim=c(60,160),on=1,col=5,lwd=1)
# add aux axis
axis(2,at=c(135,130,125,85))
axis(4,at=c(135,130,125,85))
# draw the horizontal line at the average
addSeries(xts(rep(mean(bp.xts[,1][bp.xts$High > 95]),length(index(bp.xts[bp.xts$High > 95]))),index(bp.xts[bp.xts$High > 95])),ylim=c(60,160),on=1,col=6,lwd=1)


plot(apply.weekly(apply.daily(merge(to.daily(bp.xts[bp.xts$High > 95][,1])[,2],to.daily(bp.xts[bp.xts$High > 95][,2])[,2]),mean),mean),type='p')


bp.tmp <- merge(to.daily(bp.xts[bp.xts$High > 95][,1])[,2],to.daily(bp.xts[bp.xts$High > 95][,2])[,2])
colnames(bp.tmp)[1] <- "high"
colnames(bp.tmp)[2] <- "low"
apply.weekly(bp.tmp,mean)
plot(bp.tmp ,type='p' ,ylim=c(50,170))
addSeries(xts(rep(130,length(index(bp.tmp))),index(bp.tmp)),ylim=c(50,170),on=1,col=5,lwd=1)
addSeries(xts(rep(85,length(index(bp.tmp))),index(bp.tmp)),ylim=c(50,170),on=1,col=5,lwd=1)