Manipulate timezone

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

• Set timezone to other than Tokyo, as blood pressure may be measured after the mid night.
• Setting time zone as "Bangkok" allows samples taken from 00:00 to 02:00 to be taken into the day before.

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]