2018年5月17日木曜日

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)

2018年5月8日火曜日

水上偵察機


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

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]