draw the graph to record blood pressure vol.3

The original xts package seems to have the problem to handle time based objects which belong to the timezone other than UTC. The code below to bypass the issue to normalizes TZ parameter.

1. strip index date from the original object "bp.xts" by index()
2. convert index data which contains date and time to data format by as.Date() with the parameter "tz=Sys.timezone()" or "tz=tzone(bp.xts)". Both must return "Asia/Tokyo" in my environment.
3. rejoin the numeric data and the rebuilt index by merge().

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"
apply.weekly(bp.day,mean)
# plot(apply.weekly(bp.day,mean),type='p')
#
# adjust mix-max of ylim by min() and max()
#
plot(apply.weekly(bp.day,mean),type='p',ylim=c( min(apply.weekly(bp.day,mean)[,2])-10,max(apply.weekly(bp.day,mean)[,1])+10))
#
# draw the horizontal line at 125
#
addSeries(as.xts(rep(125,length(apply.weekly(bp.day,mean)[,1])),index(apply.weekly(bp.day,mean))),on=1,ylim=c( min(apply.weekly(bp.day,mean)[,2])-10,max(apply.weekly(bp.day,mean)[,1])+10))
addSeries(as.xts(rep(85,length(apply.weekly(bp.day,mean)[,1])),index(apply.weekly(bp.day,mean))),on=1,ylim=c( min(apply.weekly(bp.day,mean)[,2])-10,max(apply.weekly(bp.day,mean)[,1])+10))
addSeries(as.xts(rep(mean(bp.xts[,2]),length(apply.weekly(bp.day,mean)[,1])),index(apply.weekly(bp.day,mean))),on=1,ylim=c( min(apply.weekly(bp.day,mean)[,2])-10,max(apply.weekly(bp.day,mean)[,1])+10),col=2)

addSeries(as.xts(rep(mean(bp.xts[,1]),length(apply.weekly(bp.day,mean)[,1])),index(apply.weekly(bp.day,mean))),on=1,ylim=c( min(apply.weekly(bp.day,mean)[,2])-10,max(apply.weekly(bp.day,mean)[,1])+10),col=2)

#
# below will use system's timezone instead of retrieving from the object.
#
bp.day <- merge(as.xts(as.vector(bp.xts[,1]),as.Date(index(bp.xts),tz=Sys.timezone())),as.vector(bp.xts[,2]))
colnames(bp.day)[1] <- "high"
colnames(bp.day)[2] <- "low"
apply.weekly(bp.day,mean)

水上偵察機

数的劣勢の下、主要仮想敵である米国海軍との決戦に挑まざるをえない帝国海軍は正面戦力の不足を挽回するため偵察能力の向上に注力した。中でも重巡洋艦に搭載された水上偵察機部隊は「艦隊の目」として期待され、最新の機材と特に選抜された優秀な水上機搭乗員が配属された。太平洋戦争の全期間を通じて大規模な海戦の有る所つねに彼ら巡洋艦水上偵察機部隊の姿があった。帝国海軍においては偵察機は単機で飛行するのが通例である。彼らは孤独のうちに任務のため飛び立ち、その多くが沈黙を守ったまま消えていったのだった。誰に見られることもなく。

How to treat residuals() output as xts object.

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

above outpus like

2000-03-01   2000-06-01   2000-09-01   2000-12-01   2001-03-01   2001-06-01   2001-09-01
-18.6280839  120.6122215  186.8549495  162.9236931   63.9004346  -15.8557805  -34.5815455
2001-12-01   2002-03-01   2002-06-01   2002-09-01   2002-12-01   2003-03-01   2003-06-01 (....)

and class are like,

> class(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] "xts" "zoo"

however, this is not able to be treated as usual xts object. then do as below.

> resi <- 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]))
> resi[,1]
                   [,1]
2000-03-01  -18.6280839
2000-06-01  120.6122215
2000-09-01  186.8549495
2000-12-01  162.9236931
2001-03-01   63.9004346
(....)

> class(resi[,1])
[1] "xts" "zoo"
>plot(resi[,1])

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

# set the periods of calculation
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.05 # 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.
# 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.
# get data from FRED and store selfregression data

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]
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]
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]


# download data from yahoo finace before this.
SP5 <- as.xts(read.zoo(read.csv("~/SP5.csv")))



### caution still in debug
if(floor(MonthsBetween(mondate(last(index(G))),mondate(Sys.Date()))) > 5){ i <- i+1}
### caution ends





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