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)

S&P500 and GDP and other macro-economical data - Case-Shiller 10-City Composite Home Price Index

getSymbols("GDP",src='FRED')
getSymbols("PAYEMS",src="FRED")
getSymbols("UNDCONTSA",src="FRED")

getSymbols('SPCS10RSA',src='FRED')
CS <- SPCS10RSA

k2k <- "2000-01-01::2017-12-01"

PAq <- apply.quarterly(PA[k2k],mean)
UCq <- apply.quarterly(UC[k2k],mean)
CSq <- apply.quarterly(CS[k2k],mean)
G <- GDP

summary(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]))

> summary(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]))

Call:

lm(formula = 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:

     Min       1Q   Median       3Q      Max 

-193.182  -33.564   -0.483   37.984  186.693 

Coefficients:

                                    Estimate Std. Error t value Pr(>|t|)    

(Intercept)                       -1.080e+03  3.794e+02  -2.846  0.00591 ** 

PAq[k2k]:UCq[k2k]                  7.692e-05  6.148e-06  12.512  < 2e-16 ***

PAq[k2k]:CSq[k2k]                 -2.019e-04  3.278e-05  -6.160 5.13e-08 ***

PAq[k2k]:UCq[k2k]:G[k2k]          -5.157e-09  4.380e-10 -11.776  < 2e-16 ***

PAq[k2k]:UCq[k2k]:CSq[k2k]        -1.639e-07  3.181e-08  -5.153 2.59e-06 ***

PAq[k2k]:G[k2k]:CSq[k2k]           2.161e-08  1.608e-09  13.442  < 2e-16 ***

PAq[k2k]:UCq[k2k]:G[k2k]:CSq[k2k]  9.053e-12  1.822e-12   4.968 5.19e-06 ***

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 75.55 on 65 degrees of freedom

Multiple R-squared:  0.9713, Adjusted R-squared:  0.9686 

F-statistic: 366.6 on 6 and 65 DF,  p-value: < 2.2e-16

> my_sp5cs(k2k,19981.13, 147810,1120,221)

[1] "m_m params! apply.quarter - UC w/ nominal GDP"

[1] 2725.425

k = [1] "2000-01-01::2017-12-01"

plot(xts(as.numeric(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]))),seq(as.Date(strsplit(k2k,"::")[[1]][1]),as.Date(strsplit(k2k,"::")[[1]][2]),by="quarters")),type='p')

Plot actual and theoretical - S&P500

plot(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]))))

plot(xts(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]))[,1],seq(as.Date(strsplit(k2k,"::")[[1]][1]),as.Date(strsplit(k2k,"::")[[1]][2]),by='quarters')),type='p')

New S&P500 = GDP and other macro-economical data

getSymbols("GDP",src='FRED')
getSymbols("PAYEMS",src="FRED")
getSymbols("UNDCONTSA",src="FRED")
getSymbols('SPCS20RSA',src='FRED')

CS <- SPCS20RSA
k2k <- "2000-01-01::2017-12-01"
PAq <- apply.quarterly(PA[kikan],mean)
UCq <- apply.quarterly(UC[kikan],mean)
CSq <- apply.quarterly(CS[kikan],mean)

summary(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:
     Min       1Q   Median       3Q      Max
-208.815  -29.865    1.545   39.597  203.971

Coefficients:
                                    Estimate Std. Error t value Pr(>|t|) 
(Intercept)                       -1.250e+03  4.070e+02  -3.071  0.00311 **
PAq[k2k]:UCq[k2k]                  8.127e-05  6.593e-06  12.327  < 2e-16 ***
PAq[k2k]:CSq[k2k]                 -2.300e-04  3.765e-05  -6.109 6.27e-08 ***
PAq[k2k]:UCq[k2k]:G[k2k]          -5.465e-09  4.813e-10 -11.355  < 2e-16 ***
PAq[k2k]:UCq[k2k]:CSq[k2k]        -1.731e-07  3.413e-08  -5.070 3.54e-06 ***
PAq[k2k]:G[k2k]:CSq[k2k]           2.503e-08  1.927e-09  12.988  < 2e-16 ***
PAq[k2k]:UCq[k2k]:G[k2k]:CSq[k2k]  9.177e-12  1.954e-12   4.697 1.41e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 80.57 on 65 degrees of freedom
Multiple R-squared:  0.9674, Adjusted R-squared:  0.9643
F-statistic:   321 on 6 and 65 DF,  p-value: < 2.2e-16



my_sp5cs(k2k,last(G[k2k]),as.vector(last(PAq[k2k])),as.vector(last(UCq[k2k])),as.vector(last(CSq[k2k])))

S&P500,GDP,PAYEMS,UNDCONTSA 2018JAN version

# set parameters.
kg <- "1992-01-01::2017-12-31"
kikan <- "1992-01-01::2017-12-31"
l <- 48  # # of months to predict
r <- 1.05 # pesumed GDP growth rate
# download data from FRE.
getSymbols("GDP",src="FRED")
getSymbols("PAYEMS",src="FRED")
getSymbols("UNDCONTSA",src="FRED")
PA <- PAYEMS
UC <- UNDCONTSA
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.
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]
my_sp5(kg,m_GDP,m_PA,m_UC)

> kikan
[1] "1992-01-01::2017-12-31"
> summary(lm(apply.quarterly(SP5[kikan],mean)[,1] ~ apply.quarterly(PA[kikan],mean) * apply.quarterly(UC[kikan],mean) * G[kikan] - apply.quarterly(UC[kikan],mean) -G[kikan] ))

Call:
lm(formula = apply.quarterly(SP5[kikan], mean)[, 1] ~ apply.quarterly(PA[kikan],
    mean) * apply.quarterly(UC[kikan], mean) * G[kikan] - apply.quarterly(UC[kikan],
    mean) - G[kikan])

Residuals:
    Min      1Q  Median      3Q     Max
-255.30  -74.08  -11.65   71.56  300.72


                                                                           Pr(>|t|) 
(Intercept)                                                                2.53e-05 ***
apply.quarterly(PA[kikan], mean)                                           1.13e-05 ***
apply.quarterly(PA[kikan], mean):apply.quarterly(UC[kikan], mean)          4.89e-10 ***
apply.quarterly(PA[kikan], mean):G[kikan]                                  3.12e-09 ***
apply.quarterly(UC[kikan], mean):G[kikan]                                   < 2e-16 ***
apply.quarterly(PA[kikan], mean):apply.quarterly(UC[kikan], mean):G[kikan]  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 117.4 on 98 degrees of freedom
Multiple R-squared:  0.9505, Adjusted R-squared:  0.9479
F-statistic: 376.1 on 5 and 98 DF,  p-value: < 2.2e-16