T test to compare before and after the new drug prescribed.

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]))
par(mfrow=c(2,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,10),col=3)
hist(as.vector(bp.bangkok["2018-06-21::"][,2]),breaks=20,xlim=c(55,100),ylim=c(0,10),col=2)
par(mfrow=c(1,1))
#
# mean value before
#
mean(last(bp.bangkok["::2018-06-20"][,2],n=length(as.vector(bp.bangkok["2018-06-21::"][,2]))))
#
# mean value after
#
mean(as.vector(bp.bangkok["2018-06-21::"][,2]))

• 2018-06-20 is the day to start to add the new drug for me.
• Comparing low blood pressure before and after the date.
• The number of sample is standardized.
• Drawing histograms to compare graphically.

.emacs

Don't forget good old days.

(global-set-key "\C-h" 'delete-backward-char)
(setenv "PATH" "/usr/bin:/bin:/usr/sbin:/sbin:/usr/local/bin")
(setq exec-path (parse-colon-path (getenv "PATH")))
; (add-to-list 'load-path "~/.emacs.d")                                       
(global-visual-line-mode)
(autoload 'wikipedia-mode "wikipedia-mode.el"
   "Major mode for editing documents in Wikipedia markup." t)
(require 'epg)
(setq epa-armor t)

; 緯度経度は小数点一桁までで指定@東京駅
(setq calendar-latitude 35.7)  ;
(setq calendar-longitude 139.8)  ;
(setq calendar-location-name "Tokyo,Japan")

H1B Visa Holder Salary

オラクル　=  Oracle America Inc

グーグル = Google Inc

で各社のH1Bビザを持ってる(??)人の給料が比較る。やはりグーグルのほうが高い。

http://h1bdata.info/index.php?em=Oracle+America+Inc&job=&city=&year=All+Years

Retrieve time data from POSIXct and filter data

bp.xts["2018-06"][strptime(format(index(bp.xts["2018-06"]),"%H:%M:%S"),"%H:%M:%S") > strptime("09:00:00","%H:%M:%S") & strptime(format(index(bp.xts["2018-06"]),"%H:%M:%S"),"%H:%M:%S") < strptime("12:00:00","%H:%M:%S")]
                    High Low Pulse
2018-06-02 10:20:00  136  91    70
2018-06-17 09:58:00  133  89    64


bp.xts[strptime(format(index(bp.xts),"%H:%M:%S"),"%H:%M:%S") > strptime("09:00:00","%H:%M:%S") & strptime(format(index(bp.xts),"%H:%M:%S"),"%H:%M:%S") < strptime("12:00:00","%H:%M:%S")]

                    High Low Pulse
2018-01-02 10:15:00  158 107    NA
2018-01-03 09:30:00  153 110    NA
2018-01-06 10:00:00  143 108    75
<skip....>
2018-05-26 09:33:00  126  85    73
2018-06-02 10:20:00  136  91    70
2018-06-17 09:58:00  133  89    64

> mean(na.omit(bp.xts[strptime(format(index(bp.xts),"%H:%M:%S"),"%H:%M:%S") > strptime("07:00:00","%H:%M:%S") & strptime(format(index(bp.xts),"%H:%M:%S"),"%H:%M:%S") < strptime("08:00:00","%H:%M:%S")])[,1])
[1] 126.0638
> mean(na.omit(bp.xts[strptime(format(index(bp.xts),"%H:%M:%S"),"%H:%M:%S") > strptime("08:00:00","%H:%M:%S") & strptime(format(index(bp.xts),"%H:%M:%S"),"%H:%M:%S") < strptime("09:00:00","%H:%M:%S")])[,1])
[1] 126.5357
> mean(na.omit(bp.xts[strptime(format(index(bp.xts),"%H:%M:%S"),"%H:%M:%S") > strptime("09:00:00","%H:%M:%S") & strptime(format(index(bp.xts),"%H:%M:%S"),"%H:%M:%S") < strptime("10:00:00","%H:%M:%S")])[,1])
[1] 132.1429
> mean(na.omit(bp.xts[strptime(format(index(bp.xts),"%H:%M:%S"),"%H:%M:%S") > strptime("10:00:00","%H:%M:%S") & strptime(format(index(bp.xts),"%H:%M:%S"),"%H:%M:%S") < strptime("11:00:00","%H:%M:%S")])[,1])
[1] 133.125

> time_base_mean <- 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]))}
> time_base_mean("08:00:00","09:00:00")
[1] 126.5357
> mapply(time_base_mean, c("08:00:00","09:00:00","10:00:00"),c("09:00:00","10:00:00","11:00:00"))
08:00:00 09:00:00 10:00:00
126.5357 132.1429 133.1250


#
# use anonymous function in line
#
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"))

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
     NaN 128.0000 125.0566 125.1176 129.8182 133.1250 132.8333 117.7000 117.4062

Upload EPS data

download eps raw data from here.

select Additional Info = index earnings

Choose necessary data and cut&paste to the new file. This is important as number only support to convert the entire original file into the separated CSV.
Save it as ".csv".

expected as below.

$nkf -Lu eps.csv
date,eps
12/31/2019,44.08
9/30/2019,42.26
<..skip...>
1988/6/30,6.22
1988/3/31,5.53


as.xts(rev(read.csv("~/eps.csv")[,2]),seq(as.Date(mondate(as.Date("1988-01-01"))),as.Date(mondate(as.Date("1988-01-01"))+ length(read.csv("~/eps.csv")[,2])*3-1),by="quarters"))

Don't forget "rev()" as eps.csv stores data in the inverse order.
Do as.Date() -> mondate() -> as.Date() conversion otherwise unable to use by="quarters"argument.

           [,1]
1988-01-01 5.53
1988-04-01 6.22
1988-07-01 6.38
1988-10-01 5.62
1989-01-01 6.74
1989-04-01 6.48
<...skip...>
            [,1]
2018-07-01 38.18
2018-10-01 40.03
2019-01-01 38.20
2019-04-01 40.22
2019-07-01 42.26
2019-10-01 44.08

calculate annual eps data with filter(). Don't forget as.vector() as as.xts() accept only vector for the 1st arg.

as.xts(as.vector(na.omit(filter(eps,rep(1,4)))),seq(as.Date(mondate(as.Date("2019-10-01"))-124*3),as.Date(mondate(as.Date("2019-10-01"))),by="quarters"))
             [,1]
1988-10-01  23.75
1989-01-01  24.96
1989-04-01  25.22
<...skip...>
2019-04-01 156.63
2019-07-01 160.71
2019-10-01 164.76

get UNDCONTSA data from www.census.gov

Download "the total number housing under construction" from here.

1. choose "New Residential Construction"
2. "1970 from 2018"
3. choose "Housing Units under Construction"
4. "Total Units" and "United States"

Open CSV file and chop unnecessary lines at header and footer. The result must be like below.

Period,Value

Jan-1970,889

Feb-1970,888

Mar-1970,890

<.........>

Mar-2018,1124

Apr-2018,1124

When file name is "~/iCloud/R/UC-1970-2018.csv" and contains data from 1970-04-01 to 2018-04-01.

as.xts(read.csv("~/iCloud/R/UC-1970-2018.csv")[,2],seq(as.Date("1970-01-01"),as.Date("2018-04-01"),by="months"))

The below calculate data length and the date of termination accordingly.

as.xts(read.csv("~/iCloud/R/UC-1970-2018.csv")[,2],as.Date(seq(mondate(as.Date("1970-01-01")),mondate(as.Date("1970-01-01"))-1+length(read.csv("~/iCloud/R/UC-1970-2018.csv")[,2]))))

           [,1]

1970-01-01  889

1970-02-01  888

1970-03-01  890

<snipet>

2018-03-01 1124

2018-04-01 1124

欧米経済の衰退と民主的世紀の終わり ―― 拡大する「権威主義的民主主義」の富とパワー

欧米経済の衰退と民主的世紀の終わり ――　拡大する「権威主義的民主主義」の富とパワー

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)

write R object into csv file.

do as below.

write.zoo(UC,file="~/ttt.csv",sep=",")

write.csv() may be useless.

write.csv(write.zoo(UC,file="~/test.csv"))

please note that below doesn't work

write.csv(write.zoo(UC),file="~/test.csv")

$HOME also makes an error

> write.csv(write.zoo(UC,file="$HOME/test.csv"))
 file(file, ifelse(append, "a", "w")) でエラー:
   コネクションを開くことができません
 追加情報:  警告メッセージ:
 file(file, ifelse(append, "a", "w")) で:
   ファイル '$HOME/test.csv' を開くことができません: No such file or directory 

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]

絵本

絵本の最後は「王子は悲しみに暮れ、国中を探しましたが、姫の姿をみたものはいませんでした」となっていているけど、あれが本当の終わりなのかな？

ちょっと続きを書いてみたよ。

「王子は王国の外、海を超え、砂漠の向こう、深い森の中をさがすため一人旅に出ることにします。

王様は王子にいいます。「王子よ。家族も国民もすべてをすてて旅に出るというのか?」
王子は答えます。「はい、私はこの世でいちばん大事なものを取り戻すため、行かなければならないのです」
王子の覚悟を知った王様は最後に声をかけます。「森のハズレの魔女を尋ねよ。手がかりが得られるかもしれぬ」

王子は森のハズレで魔女に問います。「化物に変わってしまい、行方をくらませた姫のことを知りたい」
魔女は答えます。「かわいそうな魔物の娘。お前を愛するゆえに元の姿に戻ることも出来ず化物となって世界の果てをさまようという。しかし、お前が愛する家族、お前を慕う王国の民、全てを棄てるというなら娘を救えるかもしれん。魔法の呪文を教えてやろう。娘の行き先はわたしも知らん。お前は自分の力で見つけなければならないよ」

家族のことも王国のことも忘れてしまうくらいの長い放浪の末に王子は洞窟の奥に一人の化物を見つけます。気づいた化物は王子を脅して、追い払おうとします。

この醜い姿を愛しいあの人に見られるわけにはいかない。どうか私のことは忘れて王国にかえってください。化物はそう思いました。

でも、その時王子は気づきました。化物がその懐に鏡を持っていることを。それは城で過ごした幸せな日々に贈り物として姫に与えたものでした。

『この化物が姫に違いない』

王子は森の魔女に教えられた魔法の呪文をとなえます。

その後のことは誰も知りません。漆黒の翼をもつ番いの化物が太陽を目指して高く高く飛んで行くのを見たというもの。あまりにも高く飛んだので陽の光に燃やされて灰になって地上に帰ることもなかったというもの。あるいは夕日の向こう永遠の夜を二人して今も飛び続けているというもの。


ただ一つ、王子が化物を探し当てたそのころから王国の子どもたちはオネショをしなくなったそうな。どっとはらい」