read file and normalize data to adjust TZ. for 2019

THIS ENTRY SUPERCEDES 2018-10-06.

# read data as csv format and convert to xts to plot
#
Sys.setenv(TZ=Sys.timezone())
# 
# Update below to check file.exists.
#
if(file.exists("~/~/Downloads/bp2018.csv")){
  bp2018 <- read.csv("~/Downloads/bp2018.csv")
  system("rm \"$HOME/Downloads/bp2018.csv\"")
}else{
  print("!!!FILE DOESN'T EXIST!!!!!")
}
#
# Update ends.
#
# bp2018.xts <- xts(bp2018[,c(-1,-2,-6)],as.POSIXct(paste(bp2018$Date,bp2018$Time,sep=" "),tz=Sys.timezone()),tz=Sys.timezone())
bp2018.xts <- xts(bp2018[,c(3,4,5)],as.POSIXct(paste(bp2018$Date,bp2018$Time,sep=" "),tz=Sys.timezone()),tz=Sys.timezone())
bp2019 <- read.csv("~/Downloads/bp - シート1.csv")
system("rm \"$HOME/Downloads/bp - シート1.csv\"")
bp2019.xts <- xts(bp2019[,c(3,4,5)],as.POSIXct(paste(bp2019$Date,bp2019$Time,sep=" "),tz=Sys.timezone()),tz=Sys.timezone())
bp.xts <- append(bp2018.xts,bp2019.xts)
# weekly average
apply.weekly(bp.xts[bp.xts$High > 95],mean)
#
#
# prepare data according to system timezone. "Asia/Tokyo" in most cases.
#
bp.day <- merge(as.xts(as.vector(bp.xts[,1]),as.Date(index(bp.xts),tz=tzone(bp.xts))),as.vector(bp.xts[,2]))
colnames(bp.day)[1] <- "high"
colnames(bp.day)[2] <- "low"
#
# prepare timezone 2 hours behind "Asia/Tokyo".
#
bp.bangkok <- merge(as.xts(as.vector(bp.xts[,1]),as.Date(index(bp.xts),tz="Asia/Bangkok")),as.vector(bp.xts[,2]))
colnames(bp.bangkok)[1] <- "high"
colnames(bp.bangkok)[2] <- "low"
apply.weekly(bp.bangkok,mean)


# read data as csv format and convert to xts to plot
#
Sys.setenv(TZ=Sys.timezone())
bp2018 <- read.csv("~/Downloads/bp2018.csv")
system("rm \"$HOME/Downloads/bp2018.csv\"")
# bp2018.xts <- xts(bp2018[,c(-1,-2,-6)],as.POSIXct(paste(bp2018$Date,bp2018$Time,sep=" "),tz=Sys.timezone()),tz=Sys.timezone())
bp2018.xts <- as.xts(bp2018[,c(3,4,5)],as.POSIXct(paste(bp2018$Date,bp2018$Time,sep=" "),tz=Sys.timezone()),tz=Sys.timezone())
bp2019 <- read.csv("~/Downloads/bp - シート1.csv")
system("rm \"$HOME/Downloads/bp - シート1.csv\"")
bp2019.xts <- xts(bp2019[,c(3,4,5)],as.POSIXct(paste(bp2019$Date,bp2019$Time,sep=" "),tz=Sys.timezone()),tz=Sys.timezone())
bp.xts <- append(bp2018.xts,bp2019.xts)
# weekly average
apply.weekly(bp.xts[bp.xts$High > 95],mean)
#
#
# prepare data according to system timezone. "Asia/Tokyo" in most cases.
#
bp.day <- merge(as.xts(as.vector(bp.xts[,1]),as.Date(index(bp.xts),tz=tzone(bp.xts))),as.vector(bp.xts[,2]))
colnames(bp.day)[1] <- "high"
colnames(bp.day)[2] <- "low"
#
# prepare timezone 2 hours behind "Asia/Tokyo".
#
bp.bangkok <- merge(as.xts(as.vector(bp.xts[,1]),as.Date(index(bp.xts),tz="Asia/Bangkok")),as.vector(bp.xts[,2]))
colnames(bp.bangkok)[1] <- "high"
colnames(bp.bangkok)[2] <- "low"
apply.weekly(bp.bangkok,mean)

merge all data eps spline gpuc

merge all data

result.gpuc <- 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])
result.eps <- lm(apply.quarterly(SP5[,4][k2k],mean) ~ eps_year_xts[k2k]+apply.quarterly(PA[k2k],mean)+apply.quarterly(CS[k2k],mean)+apply.quarterly(UC[k2k],mean))
SP5.result <- merge(residuals(result.gpuc),predict(result.gpuc),residuals(result.eps),predict(result.eps))

GSPC.predict <- merge(to.monthly(GSPC)[substr(k2k,11,23)],last(spline(seq(1,length(SP5.result[,1]),1),as.vector(SP5.result[,2]),n=length(SP5.result[,1])*3+1)$y,n=length(to.monthly(GSPC)[,1][substr(k2k,11,23)])),last(spline(seq(1,length(SP5.result[,1]),1),as.vector(SP5.result[,4]),n=length(SP5.result[,1])*3+1)$y,n=length(to.monthly(GSPC)[,1][substr(k2k,11,23)])),suffixes=c('','spline','eps'))


plot(merge(GSPC.predict[,4],GSPC.predict[,7],GSPC.predict[,8],GSPC.predict[,4]-GSPC.predict[,7],GSPC.predict[,4]-GSPC.predict[,8]),main="GSPC.predict[,4] vs. GSPC.predict[,7]",grid.ticks.on='months')
tmp.legend <- "Black: actual \nRed: spline\nGreen: eps"
addLegend(legend.loc = "topleft", legend.names = tmp.legend,col=3)
tmp.addTA <- as.xts(rep(2800,length(index(GSPC.predict))),index(GSPC.predict))

addSeries(tmp.addTA,on=1,col=6,lwd=1)

read file and normalize data to adjust TZ.

THIS ENTRY IS OBSOLETE. PLEASE GO TO THIS PAGE

# 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,-6)],as.POSIXct(paste(bp$Date,bp$Time,sep=" "),tz=Sys.timezone()),tz=Sys.timezone())
# weekly average
apply.weekly(bp.xts[bp.xts$High > 95],mean)
#
#
# prepare data according to system timezone. "Asia/Tokyo" in most cases.
#
bp.day <- merge(as.xts(as.vector(bp.xts[,1]),as.Date(index(bp.xts),tz=tzone(bp.xts))),as.vector(bp.xts[,2]))
colnames(bp.day)[1] <- "high"
colnames(bp.day)[2] <- "low"
#
# prepare timezone 2 hours behind "Asia/Tokyo".
#
bp.bangkok <- merge(as.xts(as.vector(bp.xts[,1]),as.Date(index(bp.xts),tz="Asia/Bangkok")),as.vector(bp.xts[,2]))
colnames(bp.bangkok)[1] <- "high"
colnames(bp.bangkok)[2] <- "low"
apply.weekly(bp.bangkok,mean)

Prepare Data --- getsymbols(), auto.arima(), as.yearmon(), as.yearqtr()

> len_mon <- 48  # # of months to predict
> r <- 1.04 # pesumed GDP growth rate
> i <- seq(2,len_mon/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.
> as.xts(forecast(auto.arima(CS["2012::"]),h=len_mon)$mean[1:len_mon],as.Date(as.yearmon(mondate(index(last(CS)))+seq(1,len_mon,1)),frac=0))[as.Date(as.yearqtr(mondate(index(last(CS)))+seq(3,len_mon,3)),frac=0)]
               [,1]
2018-07-01 225.2251
2018-10-01 227.3579
<skip>
2022-04-01 266.0765
> last(CS)
           SPCS10RSA
2018-06-01  224.8054
> as.xts(forecast(auto.arima(PA),h=len_mon)$mean[1:len_mon],as.Date(as.yearmon(mondate(index(last(PA)))+seq(1,len_mon,1)),frac=0))[as.Date(as.yearqtr(mondate(index(last(PA)))+seq(3,len_mon,3)),frac=0)]
               [,1]
2018-10-01 149684.2
2019-01-01 150292.1
<skip>
2022-07-01 158802.2
>  last(PA)
           PAYEMS
2018-08-01 149279
> last(UC)
           UNDCONTSA
2018-07-01      1122
>
as.xts(forecast(auto.arima(UC),h=len_mon)$mean[1:len_mon],as.Date(as.yearmon(mondate(index(last(UC)))+seq(1,len_mon,1)),frac=0))[as.Date(as.yearqtr(mondate(index(last(UC)))+seq(3,len_mon,3)),frac=0)]
2018-10-01 1130.871
2019-01-01 1139.742
<skip>
2022-07-01 1263.939


len_mon <- 48  # # of months to predict
gdp_g_r <- 1.04 # pesumed GDP growth rate
i <- seq(2,len_mon/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.

getSymbols("PAYEMS",src="FRED")
PA <- PAYEMS
m_PA <- as.xts(forecast(auto.arima(PA),h=len_mon)$mean[1:len_mon],as.Date(as.yearmon(mondate(index(last(PA)))+seq(1,len_mon,1)),frac=0))[as.Date(as.yearqtr(mondate(index(last(PA)))+seq(3,len_mon,3)),frac=0)]
PAq <- apply.quarterly(PA[k2k],mean)
length(PAq)

getSymbols("UNDCONTSA",src="FRED")
UC <- UNDCONTSA
m_UC <- as.xts(forecast(auto.arima(UC),h=len_mon)$mean[1:len_mon],as.Date(as.yearmon(mondate(index(last(UC)))+seq(1,len_mon,1)),frac=0))[as.Date(as.yearqtr(mondate(index(last(UC)))+seq(3,len_mon,3)),frac=0)]
UCq <- apply.quarterly(UC[k2k],mean)
length(UCq)

getSymbols('SPCS10RSA',src='FRED')
CS <- SPCS10RSA
m_CS_2012 <- as.xts(forecast(auto.arima(CS["2012::"]),h=len_mon)$mean[1:len_mon],as.Date(as.yearmon(mondate(index(last(CS)))+seq(1,len_mon,1)),frac=0))[as.Date(as.yearqtr(mondate(index(last(CS)))+seq(3,len_mon,3)),frac=0)]
CSq <- apply.quarterly(CS[k2k],mean)
length(CSq)

getSymbols("GDP",src="FRED")
G <- GDP
# m_GDP <- as.xts(as.vector(last(GDP)) * r**(i/4),d)
m_GDP <- as.xts(as.vector(last(G)) * gdp_g_r**(seq(1,len_mon/3,1)/4),as.Date(as.yearqtr(mondate(index(last(G)))+seq(3,len_mon,3)),frac=0))
kikan <- paste("1992-01-01::",as.Date(as.yearmon((mondate(index(last(G)))+2)),frac=1),sep="")
k2k <- paste("2000-01-01::",as.Date(as.yearmon((mondate(index(last(G)))+2)),frac=1),sep="")

SP5 <- as.xts(read.zoo(read.csv("~/SP5.csv")))

length(CSq)
length(UCq)
length(PAq)
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]))
my_sp5cs(k2k,m_GDP[d[1:9]],m_PA[d[1:9]],m_UC[d[1:9]],m_CS_2012[d[1:9]])
result.eps <- lm(apply.quarterly(SP5[,4][k2k],mean) ~ eps_year_xts[k2k]+apply.quarterly(PA[k2k],mean)+apply.quarterly(CS[k2k],mean)+apply.quarterly(UC[k2k],mean))
result.gpuc <- 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]))
SP5.result <- merge(residuals(result.gpuc),predict(result.gpuc),residuals(result.eps),predict(result.eps))

GSPC.predict <- merge(to.monthly(GSPC)[substr(k2k,11,23)],last(spline(seq(1,length(SP5.result[,1]),1),as.vector(SP5.result[,2]),n=length(SP5.result[,1])*3+1)$y,n=length(to.monthly(GSPC)[,1][substr(k2k,11,23)])),last(spline(seq(1,length(SP5.result[,1]),1),as.vector(SP5.result[,4]),n=length(SP5.result[,1])*3+1)$y,n=length(to.monthly(GSPC)[,1][substr(k2k,11,23)])),suffixes=c('','spline','eps'))

draw histgram -- hist function

The updated version after stopping natrix at 2018-08-09.

par(mfrow=c(4,1))
my_bp_hist_x(bp.bangkok,"::2018-06-19",2,70,0.1,0.6,20,15,55,100,4)
my_bp_hist_x(bp.bangkok,"2018-06-21::2018-07-13",2,length(bp.bangkok["2018-06-21::2018-07-13"][,2]),0.1,0.6,20,15,55,100,5)
my_bp_hist_x(bp.bangkok,"2018-07-15::2018-08-09",2,length(bp.bangkok["2018-07-15::2018-08-09"][,2]),0.1,0.6,20,15,55,100,6)
my_bp_hist_x(bp.bangkok,"2018-08-10::",2,length(bp.bangkok["2018-08-10::"][,2]),0.1,0.6,30,15,55,100,7)
par(mfrow=c(1,1))


par(mfrow=c(4,1))
my_bp_hist_x(bp.bangkok,"::2018-06-19",1,70,0.1,0.6,30,15,90,150,4)
my_bp_hist_x(bp.bangkok,"2018-06-21::2018-07-13",1,length(bp.bangkok["2018-06-21::2018-07-13"][,2]),0.1,0.6,20,15,90,150,5)
my_bp_hist_x(bp.bangkok,"2018-07-15::2018-08-09",1,length(bp.bangkok["2018-07-15::2018-08-09"][,2]),0.1,0.6,20,15,90,150,6)
my_bp_hist_x(bp.bangkok,"2018-08-10::",1,length(bp.bangkok["2018-08-10::"][,2]),0.1,0.6,30,15,90,150,7)
par(mfrow=c(1,1))

Manage Blood Pressure data

# 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,-6)],as.POSIXct(paste(bp$Date,bp$Time,sep=" "),tz=Sys.timezone()),tz=Sys.timezone())
# weekly average
apply.weekly(bp.xts[bp.xts$High > 95],mean)
#
#  draw the graph
#
plot(bp.xts[,c(1,2)][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,120,115,85))
axis(4,at=c(135,130,125,120,115,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)
addSeries(xts(rep(mean(bp.xts[,2][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)
events <- xts(c("natrix","weight","abort natrix"),as.Date(c("2018-06-20", "2018-07-14","2018-08-09")))

addEventLines(events, srt=90, pos=2,col=10)

# prepare data according to system timezone. "Asia/Tokyo" in most cases.
#
bp.day <- merge(as.xts(as.vector(bp.xts[,1]),as.Date(index(bp.xts),tz=tzone(bp.xts))),as.vector(bp.xts[,2]))
colnames(bp.day)[1] <- "high"
colnames(bp.day)[2] <- "low"
#
# prepare timezone 2 hours behind "Asia/Tokyo".
#
bp.bangkok <- merge(as.xts(as.vector(bp.xts[,1]),as.Date(index(bp.xts),tz="Asia/Bangkok")),as.vector(bp.xts[,2]))
colnames(bp.bangkok)[1] <- "high"
colnames(bp.bangkok)[2] <- "low"
apply.weekly(bp.bangkok,mean)


#
# t test between before 2018-06-20 and after 2018-06-21
#  2018-06-20 is to add natrix into prescription.
#  2018-07-14 is to start weight training.
#
t.test(as.vector(last(bp.bangkok["::2018-06-20"][,2],n=length(as.vector(bp.bangkok["2018-06-21::"][,2])))),as.vector(bp.bangkok["2018-06-21::"][,2]))
#
# t test no.2
#
t.test(as.vector(bp.bangkok["2018-06-21::2018-07-13"][,2]),as.vector(bp.bangkok["2018-07-15::"][,2]))
#
# t test no.3
#
t.test(as.vector(last(bp.bangkok["::2018-07-13"][,2],n=length(as.vector(bp.bangkok["2018-07-15::"][,2])))),as.vector(bp.bangkok["2018-07-15::"][,2]))
#
# t test no.4
#
t.test(as.vector(last(bp.bangkok["::2018-06-19"][,2],n=length(as.vector(bp.bangkok["2018-08-16::"][,2])))),as.vector(bp.bangkok["2018-08-16::"][,2]))

#
# use anonymous function in line to calculate the average of the data sorted by the time to take samples of higher pressure.
#
mapply(function(x,y){return(mean(na.omit(bp.xts[strptime(format(index(bp.xts),"%H:%M:%S"),"%H:%M:%S") > strptime(x,"%H:%M:%S") & strptime(format(index(bp.xts),"%H:%M:%S"),"%H:%M:%S") < strptime(y,"%H:%M:%S")])
[,1]))},c("05:00:00","06:00:00","07:00:00","08:00:00","09:00:00","10:00:00","11:00:00","23:00:00","00:00:00"),c("05:00:00","06:59:00","07:59:00","08:59:00","09:59:00","10:59:00","11:59:00","23:59:00","00:59:00"))
#
# mean value of the lower prssure.
#
mapply(function(x,y){return(mean(na.omit(bp.xts[strptime(format(index(bp.xts),"%H:%M:%S"),"%H:%M:%S") > strptime(x,"%H:%M:%S") & strptime(format(index(bp.xts),"%H:%M:%S"),"%H:%M:%S") < strptime(y,"%H:%M:%S")])
[,2]))},c("05:00:00","06:00:00","07:00:00","08:00:00","09:00:00","10:00:00","11:00:00","23:00:00","00:00:00"),c("05:00:00","06:59:00","07:59:00","08:59:00","09:59:00","10:59:00","11:59:00","23:59:00","00:59:00"))
#
#  count the number of samples
#
mapply(function(x,y){return(length(na.omit(bp.xts[strptime(format(index(bp.xts),"%H:%M:%S"),"%H:%M:%S") > strptime(x,"%H:%M:%S") & strptime(format(index(bp.xts),"%H:%M:%S"),"%H:%M:%S") < strptime(y,"%H:%M:%S")])
[,1]))},c("05:00:00","06:00:00","07:00:00","08:00:00","09:00:00","10:00:00","11:00:00","23:00:00","00:00:00"),c("05:00:00","06:59:00","07:59:00","08:59:00","09:59:00","10:59:00","11:59:00","23:59:00","00:59:00"))


par(mfrow=c(3,1))
hist(as.vector(last(bp.bangkok["::2018-06-20"][,2],n=length(as.vector(bp.bangkok["2018-06-21::"][,2])))),breaks=20,xlim=c(55,100),ylim=c(0,15),col=3)
# axis(side=2, pos=84,labels=F) 
axis(side=2, pos=round(mean(bp.bangkok["::2018-06-20"][,2])),labels=F) 
hist(as.vector(bp.bangkok["2018-06-21::2018-07-13"][,2]),breaks=20,xlim=c(55,100),ylim=c(0,15),col=2)
axis(side=2, pos=round(mean(bp.bangkok["2018-06-21::2018-07-13"][,2])),labels=F) 
hist(as.vector(bp.bangkok["2018-07-15::"][,2]),breaks=10,xlim=c(55,100),ylim=c(0,15),col=4)
axis(side=2, pos=round(mean(bp.bangkok["2018-07-15::"][,2])),labels=F) 
graph_dim <- par('usr') 
text( (graph_dim[1] + graph_dim[2]) / 2   ,(graph_dim[4] - graph_dim[3]) * 0.70 + graph_dim[3] ,paste("mean",mean(bp.bangkok["2018-07-15::"][,2]),sep="="),adj=c(0,0))
text( (graph_dim[1] + graph_dim[2]) / 2   ,(graph_dim[4] - graph_dim[3]) * 0.75 + graph_dim[3] ,paste("sd",sd(bp.bangkok["2018-07-15::"][,2]),sep="="),adj=c(0,0))
par(mfrow=c(1,1))

Prepare Data and set parameters.

kg <- "1992-01-01::2018-03-31"
kikan <- "1992-01-01::2018-03-31"
k2k <- "2000-01-01::2018-03-31"   #GPUC model is based data only from 2000 because of cs index availability
# prepare parameters.
l <- 48  # # of months to predict
r <- 1.04 # pesumed GDP growth rate
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.

getSymbols("GDP",src="FRED")
G <- GDP
m_GDP <- as.xts(as.vector(last(GDP)) * r**(i/4),d)

getSymbols("PAYEMS",src="FRED")
PA <- PAYEMS
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]
PAq <- apply.quarterly(PA[k2k],mean)

getSymbols('SPCS10RSA',src='FRED')
CS <- SPCS10RSA
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]
CSq <- apply.quarterly(CS[k2k],mean)

getSymbols("UNDCONTSA",src="FRED")
# UC <- UNDCONTSA
# comment out here until UNDCONTSA update is resumed in FRED.
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]
UCq <- apply.quarterly(UC[k2k],mean)

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]

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')

Calculate NIKKEI225 vol.3

k3 <-  "2007-01-01::2018-02-07"

# or

getSymbols("^GSPC")
k3 <- paste("2007-01-01", index(last(GSPC)),sep="::")
k3

# download other data
getSymbols("NIKKEI225",src="FRED") # download nikkei 225
# getSymbols("DEXJPUS", src = "FRED")
# getSymbols("YJUSDJPY",src="yahooj")
getSymbols('YJUSDJPY', src="yahooj",auto.assign=TRUE)
N225 <- NIKKEI225
# result <- summary(lm(to.monthly(N225[k3])[,4] ~  to.monthly(GSPC[k3])[,4] + to.monthly(DEXJPUS[k3])[,4]))  
result_nikkei <- lm(to.monthly(N225[k3])[,4] ~  to.monthly(GSPC[k3])[,4] + to.monthly(YJUSDJPY[k3])[,4])
result_nikkei$coefficients[2]*last(GSPC)[,4]+result_nikkei$coefficients[3]*as.vector(last(YJUSDJPY)[,4])+result_nikkei$coefficients[1]
# plot(merge(as.xts(predict(result_nikkei),index(residuals(result_nikkei))),to.monthly(N225[k3])[,4],residuals(result_nikkei)))
beep(2)
plot(merge(as.xts(predict(result_nikkei),index(residuals(result_nikkei))),to.monthly(N225[k3])[,4],residuals(result_nikkei)),main=paste(paste("NIKKEI225 =",round(result_nikkei$coefficients[2],4),"* GSPC +",round(result_nikkei$coefficients[3],2),"*USDJPY +",round(result_nikkei$coefficients[1],2))))

tmp.legend <- paste("R Squared is ",round(summary(result_nikkei)$r.squared,4)," \n","DF is ",round(summary(result_nikkei)$df[2],0),sep=' ')

addLegend(legend.loc = "topleft", legend.names = tmp.legend,col=3)

> result

Call:
lm(formula = to.monthly(N225[k3])[, 4] ~ to.monthly(GSPC[k3])[,
    4] + to.monthly(YJUSDJPY[k3])[, 4])

Residuals:
     Min       1Q   Median       3Q      Max
-1875.41  -473.74   -12.89   507.10  1855.32

Coefficients:
                                Estimate Std. Error t value Pr(>|t|) 
(Intercept)                   -9864.4555   499.7539  -19.74   <2e-16 ***
to.monthly(GSPC[k3])[, 4]         4.8991     0.1804   27.16   <2e-16 ***
to.monthly(YJUSDJPY[k3])[, 4]   158.2244     6.1444   25.75   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 792.4 on 131 degrees of freedom
Multiple R-squared:  0.9655, Adjusted R-squared:  0.965
F-statistic:  1832 on 2 and 131 DF,  p-value: < 2.2e-16


summary(lm(apply.weekly(N225[k3],mean) ~  apply.weekly(GSPC[k3][,4],mean) + apply.weekly(YJUSDJPY[k3][,4],mean)))   

> summary(lm(apply.weekly(N225[k3],mean) ~  apply.weekly(GSPC[k3][,4],mean) + apply.weekly(YJUSDJPY[k3][,4],mean)))   

Call:
lm(formula = apply.weekly(N225[k3], mean) ~ apply.weekly(GSPC[k3][,
    4], mean) + apply.weekly(YJUSDJPY[k3][, 4], mean))

Residuals:
     Min       1Q   Median       3Q      Max
-1992.66  -510.50    51.79   469.11  2524.40

Coefficients:
                                        Estimate Std. Error t value Pr(>|t|) 
(Intercept)                           -1.017e+04  2.593e+02  -39.24   <2e-16 ***
apply.weekly(GSPC[k3][, 4], mean)      4.800e+00  9.568e-02   50.17   <2e-16 ***
apply.weekly(YJUSDJPY[k3][, 4], mean)  1.626e+02  3.180e+00   51.11   <2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 754.9 on 445 degrees of freedom
  (132 observations deleted due to missingness)
Multiple R-squared:  0.9674, Adjusted R-squared:  0.9673
F-statistic:  6612 on 2 and 445 DF,  p-value: < 2.2e-16

# plot(merge(na.omit(apply.weekly(N225[k3],mean))[2:length(na.omit(apply.weekly(N225[k3],mean)))],predict(lm(apply.weekly(N225[k3],mean) ~  apply.weekly(GSPC[k3][,4],mean) + apply.weekly(YJUSDJPY[k3][,4],mean)))))

result <- lm(apply.weekly(N225[k3],mean) ~  apply.weekly(GSPC[k3][,4],mean) + apply.weekly(YJUSDJPY[k3][,4],mean)) 
tmp.main <- paste("NIKKEI225 =",round(result$coefficients[2],4),"* GSPC +",round(result$coefficients[3],2),"*USDJPY +",round(result$coefficients[1],2))
nikkei_week <- na.omit(apply.weekly(N225[k3],mean))
predict_nikkei_week <- predict(lm(apply.weekly(N225[k3],mean) ~  apply.weekly(GSPC[k3][,4],mean) + apply.weekly(YJUSDJPY[k3][,4],mean)))
if(length(nikkei_week) > length(predict_nikkei_week)){
  plot(merge(last(nikkei_week,n=length(predict_nikkei_week),predict_nikkei_week)),main=tmp.main)
}else{
  plot(merge(nikkei_week,last(predict_nikkei_week,n=length(nikkei_week))),main=tmp.main) 
}
# add legend at topleft
#
#
addLegend(legend.loc = "topleft", legend.names = paste("R Squared",round(summary(result)$r.squared,4),sep=':'),col=2)

or

tmp.legend <- paste("R Squared is ",round(summary(result)$r.squared,4)," \n","DF is ",round(summary(result)$df[2],0),sep=' ')
addLegend(legend.loc = "topleft", legend.names = tmp.legend,col=3)

In spite of the recent steep declines, Nikkei225 is still in the premium side.

use grid.ticks.on='months' to show grid on monthly basis.