2020年12月25日金曜日

R 4.0 ggplot

scale_fill_hue(name='regions',labels= as.character(unique(w[,5])) )

ラベルの使用順序はscale で明示すること。これがないとalphabeticalにソートされてしまう。

上記だけをやると列名とラベルが不整合を起こすので、列名にも数字を先頭に付加してソート順序を制御しなくてはならない。


  for(i in seq(1,47,1)){ colnames(mdf)[i] <-  (paste(sprintf("%02d",i),colnames(mdf)[i],sep=""))}


reshape2はもう使えない。したがってmelt() も使えない。

melt(data=mdf, id.vars="t", measure.vars=as.character(unique(w[,5]))) は

  df.melt <- mdf  %>% tidyr::gather(variable,value,as.character(colnames(mdf)[-48]))

になる。

以下をインストールしておくこと。

  • tidyr                                Tidy Messy Data
  • dplyr                                A Grammar of Data Manipulation


unloadするときは。

detach("package:dplyr", unload=TRUE)
detach("package:tidyr", unload=TRUE)



2020年11月14日土曜日

EPS 2020NOV14

 > eps_year_xts["2020::"]
             [,1]
2020-01-01 116.33
2020-04-01  99.23
2020-07-01  93.03
2020-10-01  87.16
2021-01-01 101.75
2021-04-01 119.53
2021-07-01 128.19
2021-10-01 135.75


$ tac eps.txt | awk '{gsub("\\$","",$NF);print "eps_year_xts[\"2019::\"]["NR"] <- "$NF}'

eps_year_xts["2019::"][1] <- 134.39
eps_year_xts["2019::"][2] <- 135.27
eps_year_xts["2019::"][3] <- 132.90
eps_year_xts["2019::"][4] <- 139.47
eps_year_xts["2019::"][5] <- 116.33
eps_year_xts["2019::"][6] <- 99.23
eps_year_xts["2019::"][7] <- 96.47
eps_year_xts["2019::"][8] <- 92.18
eps_year_xts["2019::"][9] <- 106.60
eps_year_xts["2019::"][10] <- 125.68
eps_year_xts["2019::"][11] <- 132.68
eps_year_xts["2019::"][12] <- 140.24

2020年10月25日日曜日

EPS 2020OCT25

 


>   eps_year_xts["2020::"]
             [,1]
2020-01-01 116.33
2020-04-01  99.16
2020-07-01  93.44
2020-10-01  89.68
2021-01-01 110.65
2021-04-01 129.66
2021-07-01 139.84
2021-10-01 148.29

$   tac eps.txt | awk '{gsub("\\$","",$NF);print "eps_year_xts[\"2019::\"]["NR"] <- "$NF}'

eps_year_xts["2019::"][1] <- 134.39
eps_year_xts["2019::"][2] <- 135.27
eps_year_xts["2019::"][3] <- 132.90
eps_year_xts["2019::"][4] <- 139.47
eps_year_xts["2019::"][5] <- 116.33
eps_year_xts["2019::"][6] <- 99.23
eps_year_xts["2019::"][7] <- 93.03
eps_year_xts["2019::"][8] <- 87.16
eps_year_xts["2019::"][9] <- 101.75
eps_year_xts["2019::"][10] <- 119.53
eps_year_xts["2019::"][11] <- 128.19
eps_year_xts["2019::"][12] <- 135.75


for reference

$    cat eps.txt 

12/31/2021 $44.69 $37.23 21.02 25.44 $164.30 $135.75
9/30/2021 $42.32 $36.45 22.20 26.94 $155.59 $128.19
6/30/2021 $40.10 $35.60 23.54 28.89 $146.68 $119.53
3/31/2021 $37.19 $26.47 25.89 33.94 $133.37 $101.75
12/31/2020 $35.98 $29.66 29.85 39.62 $115.68 $87.16
9/30/2020 (25.0%) 3363.00 $33.41 $27.79 29.05 37.12 $118.88 $93.03
6/30/2020 3100.29 $26.79 $17.83 24.75 31.24 $125.28 $99.23
3/31/2020 2584.59 $19.50 $11.88 18.64 22.22 $138.63 $116.33
12/31/2019 3230.78 $39.18 $35.53 20.56 23.16 $157.12 $139.47
9/30/2019  2976.74 $39.81 $33.99 19.46 22.40 $152.97 $132.90
6/30/2019 2941.76 $40.14 $34.93 19.04 21.75 $154.54 $135.27
3/31/2019 2834.40 $37.99 $35.02 18.52 21.09 $153.05 $134.39

2020年10月9日金曜日

awk tr sed mondate CLI download

1. download cli from this url 

1.1 change downloaded file name to CLI3.csv

1.2 ~/R/R2/index/cli_download.sed  run at the terminal.


sed -n '/USA/p' CLI3.csv |awk -F, '{print $6"-01,"$7}'  |sed 's/\"//g' |awk 'BEGIN{print "DATE,DATA"}{print $0}' > usa.csv

# extract OECD entries and exclude OECDE

sed -n '/OECD[^E]/p' CLI3.csv |awk -F, '{print $6"-01,"$7}'  |sed 's/\"//g' |awk 'BEGIN{print "DATE,DATA"}{print $0}' > oecd.csv

# 

sed -n '/CHN/p' CLI3.csv |awk -F, '{print $6"-01,"$7}'  |sed 's/\"//g' |awk 'BEGIN{print "DATE,DATA"}{print $0}' > chn.csv

# 

sed -n '/EA19/p' CLI3.csv |awk -F, '{print $6"-01,"$7}'  |sed 's/\"//g' |awk 'BEGIN{print "DATE,DATA"}{print $0}' > ea19.csv

#

sed -n '/JPN/p' CLI3.csv |awk -F, '{print $6"-01,"$7}'  |sed 's/\"//g' |awk 'BEGIN{print "DATE,DATA"}{print $0}' > jpn.csv



1.3 run cli_load.r in R.


2.when cli is downloaded form this url 


2.1. do as below in terminal change FILE_NAME accordingly.

grep "OECD - Total" <FILE_NAME>.csv  | grep LOLITOAA |gawk -F, 'BEGIN{print "c("}{print $(NF-2)","}END{print ")  " }' | tr -d "\n" | sed 's/,)/)/'


output is like 


c(100.1448,100.0088,99.88226,99.7745,99.68903,99.62392,99.56622,99.5041,99.44289,99.39536,99.37106,99.38042,99.41794,99.46317,99.49193,99.47968,99.41767,97.69222,93.17708,94.81512,97.03869,98.08919,98.53291,98.79954) 


> length(c(100.1448,100.0088,99.88226,99.7745,99.68903,99.62392,99.56622,99.5041,99.44289,99.39536,99.37106,99.38042,99.41794,99.46317,99.49193,99.47968,99.41767,97.69222,93.17708,94.81512,97.03869,98.08919,98.53291,98.79954) )

[1] 24


# start and end dates are fixed. caution!

seq( as.Date(mondate(Sys.Date())-24)-day(mondate(Sys.Date()))+1,as.Date(mondate(Sys.Date())-1),by='months' )


>seq( as.Date(mondate(Sys.Date())-24)-day(mondate(Sys.Date()))+1,as.Date(mondate(Sys.Date())-1),by='months' )

 [1] "2018-10-01" "2018-11-01" "2018-12-01" "2019-01-01" "2019-02-01" "2019-03-01" "2019-04-01" "2019-05-01" "2019-06-01"

[10] "2019-07-01" "2019-08-01" "2019-09-01" "2019-10-01" "2019-11-01" "2019-12-01" "2020-01-01" "2020-02-01" "2020-03-01"

[19] "2020-04-01" "2020-05-01" "2020-06-01" "2020-07-01" "2020-08-01" "2020-09-01"

 

2.2. do as below in R

as.xts(c(100.1448,100.0088,99.88226,99.7745,99.68903,99.62392,99.56622,99.5041,99.44289,99.39536,99.37106,99.38042,99.41794,99.46317,99.49193,99.47968,99.41767,97.69222,93.17708,94.81512,97.03869,98.08919,98.53291,98.79954), seq( as.Date(mondate(Sys.Date())-24)-day(mondate(Sys.Date()))+1,as.Date(mondate(Sys.Date())-1),by='months' ))


2.3. alternatively do as below 

grep "OECD - Total" <FILE_NAME>.csv  | grep LOLITOAA |gawk -F, 'BEGIN{print "as.xts(c("}{print $(NF-2)","}END{print "), seq( as.Date(mondate(Sys.Date())-24)-day(mondate(Sys.Date()))+1,as.Date(mondate(Sys.Date())-1),by=\"months\" ))  " }' | tr -d "\n" | sed 's/,)/)/'



2020年10月4日日曜日

mac addr DHCP LAN

 


  • 58:E2:8F:22:5F:C0       iphone se
  • F8:4E:73:02:B7:6C      (iPhone11onomoto) via br0 ?
  • 34:A8:EB:81:58:4A       ipad mini
  • F8:FF:C2:65:B5:0C       macbook pro 16inch
  • 78:4F:43:98:2F:8C       macbook pro 15inch
  • 00:1E:C2:A7:D3:EB      macbook kuro
  • 18:C2:BF:16:51:90      LS5110 NAS  wired LAN
  • 82:BA:85:04:AB:ED      iphone 11 pro
  • 74:75:48:1E:CD:4F       kindle  from mac address
  • 00:F7:6F:CF:4E:E6 airmac express
  • 58:55:CA:30:E6:4C      apple tv
  • 50:c4:dd:e9:70:b0        router

身元不明macアドレス
  • 192.168.11.102        74:75:48:1E:CD:4F       31:15:20        kindle 
  • 192.168.11.110 72:8C:56:0F:7B:73

2020年10月3日土曜日

bash for and while

 

全てのcsvを0 length にする。

for i in *.csv ; do  cat /dev/null >  $i; done


無限ループ。「while : (コロン)」で無限ループになる。

$  while : ; do echo 'hello'; sleep 10 ; done

2020年9月1日火曜日

行列の計算 matrix()

 

> 1:15
 [1]  1  2  3  4  5  6  7  8  9 10 11 12 13 14 15

> matrix(1:15,ncol=3)
     [,1] [,2] [,3]
[1,]    1    6   11
[2,]    2    7   12
[3,]    3    8   13
[4,]    4    9   14
[5,]    5   10   15

> diff(matrix(1:15,ncol=3))
     [,1] [,2] [,3]
[1,]    1    1    1
[2,]    1    1    1
[3,]    1    1    1
[4,]    1    1    1


> c(1:7,NA,9:15)
 [1]  1  2  3  4  5  6  7 NA  9 10 11 12 13 14 15
> w <- matrix(c(1:7,NA,9:15), ncol=5)
> w
     [,1] [,2] [,3] [,4] [,5]
[1,]    1    4    7   10   13
[2,]    2    5   NA   11   14
[3,]    3    6    9   12   15

> v <- as.vector(w)
> v[is.na(v)]
[1] NA
> v[is.na(v)] <- 8
> matrix(v,ncol=dim(w)[2])
     [,1] [,2] [,3] [,4] [,5]
[1,]    1    4    7   10   13
[2,]    2    5    8   11   14
[3,]    3    6    9   12   15
> w <- matrix(v,ncol=dim(w)[2])
> w
     [,1] [,2] [,3] [,4] [,5]
[1,]    1    4    7   10   13
[2,]    2    5    8   11   14
[3,]    3    6    9   12   15

EPS 2020SEP01

 


> eps_year_xts["2020::"]
             [,1]
2020-01-01 116.33
2020-04-01 100.29
2020-07-01  94.83
2020-10-01  92.54
2021-01-01 112.31
2021-04-01 128.51
2021-07-01 138.90
2021-10-01 146.96


$  tac eps.txt | awk '{gsub("\\$","",$NF);print "eps_year_xts[\"2019::\"]["NR"] <- "$NF}'
eps_year_xts["2019::"][1] <- 134.39
eps_year_xts["2019::"][2] <- 135.27
eps_year_xts["2019::"][3] <- 132.90
eps_year_xts["2019::"][4] <- 139.47
eps_year_xts["2019::"][5] <- 116.33
eps_year_xts["2019::"][6] <- 99.16
eps_year_xts["2019::"][7] <- 93.44
eps_year_xts["2019::"][8] <- 89.68
eps_year_xts["2019::"][9] <- 110.65
eps_year_xts["2019::"][10] <- 129.66
eps_year_xts["2019::"][11] <- 139.84
eps_year_xts["2019::"][12] <- 148.29

2020年8月16日日曜日

plot3d color col=

 


Don't forget "col="  in the parameters. it's not scatterplot3d!!


library(rgl)

len <- dim(mdf)[1]

x3d <- c()

for(i in seq(1,47,1)){ x3d <- append(x3d,rep(i,len))}

z3d <- c()

for(i in seq(1,47,1)){ z3d <- append(z3d,dmdf[,i])}

y3d <- as.integer(gsub('-','',as.character(dmdf$t)))

# aichi 23, osaka 27, fukuoka 40, tokyo 13.

# color <- c(rep("grey",12*len),rep("orange",len),rep("grey",9*len),rep("pink",len),rep("grey",3*len),rep("yellow",len),rep("grey",12*len),rep("green",len),rep("grey",6*len),rep("red",len))

par(bg = 'grey25', fg = 'white',col.axis = 'white',col.lab='white')

y3d <- seq(1,len,1)

# plot3d(x3d,rep(y3d,47),z3d,col = color, type = "h", pch = " ",zlim=c(0, max(dmdf[,-48])))

plot3d(x3d,rep(y3d,47),z3d,col = as.vector(matrix(rep(rainbow(47),len),ncol=47,byrow=t)), type = "h", pch = " ",zlim=c(0, max(dmdf[,-48])))




==


2020年8月11日火曜日

scatterplot3d

 use same data.frame mdf as the previous entry.


library("scatterplot3d")

len <- dim(mdf)[1]

x3d <- c()

for(i in seq(1,47,1)){ x3d <- append(x3d,rep(i,len))}

z3d <- c()

for(i in seq(1,47,1)){ z3d <- append(z3d,mdf[,i])}

# scatterplot3d(x3d,rep(mdf$t,47),z3d,highlight.3d = TRUE, type = "h", pch = " ",zlim=c(0,500))

y3d <- as.integer(gsub('-','',as.character(mdf$t)))

# scatterplot3d(x3d,rep(y3d,47),z3d,highlight.3d = TRUE, type = "h", pch = " ",zlim=c(0,500))

color <- c(rep("black",12*len),rep("orange",len),rep("black",13*len),rep("yellow",len),rep("black",12*len),rep("green",len),rep("black",7*len))

# scatterplot3d(x3d,rep(y3d,47),z3d,highlight.3d = TRUE, type = "h", pch = " ",zlim=c(0,30),angle = 65)

scatterplot3d(x3d,rep(y3d,47),z3d,color, type = "h", pch = " ",zlim=c(0, max(mdf[,-48])),angle = 65,col.grid="grey")



color <- c(rep("grey",12*len),rep("orange",len),rep("grey",13*len),rep("yellow",len),rep("grey",12*len),rep("green",len),rep("grey",7*len))
par(bg = 'grey25', fg = 'white',col.axis = 'white',col.lab='white')
scatterplot3d(x3d,rep(y3d,47),z3d,color, type = "h", pch = " ",zlim=c(0, max(mdf[,-48])),angle = 65)


2020年8月10日月曜日

scatterplot3d

when mdf is like below

head(mdf,2)

  Hokkaido Aomori Iwate Miyagi Akita Yamagata Fukushima Ibaraki Tochigi Gunma Saitama Chiba Tokyo Kanagawa Niigata Toyama

1       10      0     0      0     0        0         0       0       0     0       6     2     2        0       3      0

2        8      0     0      0     0        0         0       0       0     2       3     2     2        3       0      0

  Ishikawa Fukui Yamanashi Nagano Gifu Shizuoka Aichi Mie Shiga Kyoto Osaka Hyogo Nara Wakayama Tottori Shimane Okayama

1        0     0         0      0    0        1     7   5     0     2     9     8    0        0       0       0       0

2        0     0         0      0    0        0     3   0     0     0     3    13    0        0       0       0       0

  Hiroshima Yamaguchi Tokushima Kagawa Ehime Kochi Fukuoka Saga Nagasaki Kumamoto Oita Miyazaki Kagoshima Okinawa

1         0         0         0      0     0     0       0    0        0        0    0        0         0       0

2         0         0         0      0     0     0       0    1        0        0    0        0         0       0

           t

1 2020-03-12

2 2020-03-13


do as below

library("scatterplot3d")

len <- dim(mdf)[1]

x3d <- c()

for(i in seq(1,47,1)){ x3d <- append(x3d,rep(i,len))}

z3d <- c()

for(i in seq(1,47,1)){ z3d <- append(z3d,mdf[,i])}

# scatterplot3d(x3d,rep(mdf$t,47),z3d,highlight.3d = TRUE, type = "h", pch = " ",zlim=c(0,500))

y3d <- as.integer(gsub('-','',as.character(mdf$t)))

# scatterplot3d(x3d,rep(y3d,47),z3d,highlight.3d = TRUE, type = "h", pch = " ",zlim=c(0,500))

scatterplot3d(x3d,rep(y3d,47),z3d,highlight.3d = TRUE, type = "h", pch = " ",zlim=c(0,500),angle = 65)





2020年7月28日火曜日

R 無限ループ



library(reshape2)
library(ggplot2)
library(audio)
for(i in seq(1,36,1)){
  source("../../Dropbox/R-script/covid/tokyo_effective_repro.r")
  # source("../../Dropbox/R-script/covid/em_region.r")
  # source("../../Dropbox/R-script/covid/tokyo_age_split.r")
  source("../../Dropbox/R-script/covid/plot_effective_repro.r")
  # source("../../Dropbox/R-script/covid/tokyo_vs_others.r")
  # source("../../Dropbox/R-script/covid/em_region.r")

  print(Sys.time())
{if(sum(mdf[(length(rownames(mdf))-6):length(rownames(mdf)),i])>1){cat(colnames(mdf)[i]);cat(" ");cat(sum(mdf[(length(rownames(mdf))-6):length(rownames(mdf)),i]));t <- t+sum(mdf[(length(rownames(mdf))-6):length(rownames(mdf)),i]) ;cat("\n")}};cat(paste("total ",t,"\n",sep=''))

  wait(3600)
}

2020年7月20日月曜日

EPS 2020JUL20

> eps_year_xts["2020::"]
             [,1]
2020-01-01 116.43
2020-04-01 101.44
2020-07-01  94.97
2020-10-01  92.06
2021-01-01 113.24
2021-04-01 127.39
2021-07-01 137.85
2021-10-01 145.28


delete 2 lines in the middle of eps.txt which include 'ACTUAL' and 'price'.


~$ tac eps.txt | awk '{gsub("\\$","",$NF);print "eps_year_xts[\"2019::\"]["NR"] <- "$NF}'
eps_year_xts["2019::"][1] <- 134.39
eps_year_xts["2019::"][2] <- 135.27
eps_year_xts["2019::"][3] <- 132.90
eps_year_xts["2019::"][4] <- 139.47
eps_year_xts["2019::"][5] <- 116.33
eps_year_xts["2019::"][6] <- 100.29
eps_year_xts["2019::"][7] <- 94.83
eps_year_xts["2019::"][8] <- 92.54
eps_year_xts["2019::"][9] <- 112.31
eps_year_xts["2019::"][10] <- 128.51
eps_year_xts["2019::"][11] <- 138.90
eps_year_xts["2019::"][12] <- 146.96

2020年7月14日火曜日

年齢別グラフ g <- ggplot(NULL)



  • 東京の感染者数ヒストリカルデータから日次の年齢別ヒストグラムを作る
  • その上に重症化可能性指数を線グラフで重ねる。
  • g <- ggplot(NULL)を使用してベースオブジェクトを生成し、それに対しgeom_bar(),geom_line()を行ってグラフをオーバーレイ表示する。
  • 重症化可能性指数は大体最終の死亡者の10倍になるように若い人は人数に0.02
    50代は 0.04 60代は 0.15, 70代は 0.56 それ以上は1 の係数をかけて加算
  • 事前に最新のCSV をダウンロードしておくこと。
    • curl <- "https://stopcovid19.metro.tokyo.lg.jp/data/130001_tokyo_covid19_patients.csv"
    • cdestfile <- "~/R/R2/covid/tokyo.csv"
    • download.file(curl,cdestfile)


func <- function(x1,x2,x3,x4,x5,x6,x7,x8){
  return(x1*0.02+x2*0.02+x3*0.02+x4*0.02+x5*0.04+x6*0.15+x7*0.56+x8*1)
}

length_graph <- length(seq(as.Date("2020-03-20"),Sys.Date(),by='days'))
w <- read.csv("~/R/R2/covid/tokyo.csv")
y <- as.xts(as.numeric(substr(w[,9],1,2)),as.Date(w[,5]))
# apply.daily(as.xts(rep(1,length(y[y[,1] == 10])),as.Date(index(y[y[,1] == 10]))),sum)

# v <- c()
# seq(as.Date(w[1,5]),Sys.Date(),by='days')
v <- as.xts(rep(1,length( seq(as.Date(w[1,5]),Sys.Date(),by='days') )), seq(as.Date(w[1,5]),Sys.Date(),by='days'))
for(i in seq(10,80,10)){
  if(i < 80){
    v<- merge(v,  apply.daily(as.xts(rep(1,length(y[y[,1] == i])),as.Date(index(y[y[,1] == i]))),sum))
  }else{  # for the case the sample is more than 80 yrs old.
    v<- merge(v,  apply.daily(as.xts(rep(1,length(y[y[,1] >= i])),as.Date(index(y[y[,1] >= i]))),sum))
  }
}
# v<- merge(v,  apply.daily(as.xts(rep(1,length(y[y[,1] > 80 ])),as.Date(index(y[y[,1] >80]))),sum))
v <- v[,-1]
for(i in seq(1,8,1)){
  colnames(v)[i] <- as.character(i*10)
}

# replace NA with ZERO.
for( i in seq(1,length(colnames(v)),1)) {
  for(j in seq(1,length(index(v)),1)) {
      if(is.na(v[j,i])){  v[j,i] <- 0}
  }
  # print(v[,i])
}

# for( i in length(colnames(v))){
#  replace(v[,i], which(is.na(v[,i])), 0)
# }

# apply(v[,c(6,7,8)],1,sum)
w <- merge(as.xts(apply(v[,c(1,2,3,4,5)],1,sum) ,index(v)),as.vector(apply(v[,c(6,7,8)],1,sum) ))

colnames(w)[1] <- "lessthan60"
colnames(w)[2] <- "over60"

df <- data.frame(t=last(index(w),length_graph),
                 # o=last(diff(all[,1]-all[,3],t=all[,3]),length_graph),
                 o=last(w[,1],length_graph),
                 k=last(w[,2],length_graph)
                 # or=last(all[,2],length_graph)*multi,
                 # kr=last(all[,2],length_graph)*multi
               )
# df.melt <- melt(data=df, id.vars="t", measure.vars=c("lessthan60", "over60"))


df.melt <- melt(data=df, id.vars="t", measure.vars=c(colnames(w)[1],colnames(w)[2]))
# head(df.melt)
df <- df.melt
               # g <- ggplot(x, aes(x = t, y = value, fill = variable))
               # # g <- ggplot(x, aes(x = t, y = d))
               # g <- g + geom_bar(stat = "identity")
               # # g <- g + scale_fill_nejm()
g <- ggplot(df, aes(x = t, y = value, fill = variable))
g <- g + geom_bar(stat = "identity")
# plot(g)

df <- data.frame(t=last(index(v),length_graph),
                 # o=last(diff(all[,1]-all[,3],t=all[,3]),length_graph),
                 a=last(v[,1],length_graph),
                 b=last(v[,2],length_graph),
                 c=last(v[,3],length_graph),
                 d=last(v[,4],length_graph),
                 e=last(v[,5],length_graph),
                 f=last(v[,6],length_graph),
                 g=last(v[,7],length_graph),
                 h=last(v[,8],length_graph)
                 # or=last(all[,2],length_graph)*multi,
                 # kr=last(all[,2],length_graph)*multi
               )

for(i in seq(1,8,1)){
                 colnames(df)[i+1] <- as.character(i*10)
}

# df.melt <- melt(data=df, id.vars="t", measure.vars=c("X10", "X20", "X30", "X40", "X50", "X60", "X70", "X80"))
df.melt <- melt(data=df, id.vars="t", measure.vars=c("10", "20", "30", "40", "50", "60", "70", "80"))
# head(df.melt)
df <- df.melt
# in order to overlayer graph use ggplot(NULL) to create base object.
g <- ggplot(NULL)
# g <- ggplot(df, aes(x = t, y = value, fill = variable))
g <- g + scale_fill_brewer(palette="Spectral",na.value = "black",name = "age group", direction=-1,labels = c("=<19",">=20",">=30",">=40",">=50",">=60",">=70",">=80"))
g <- g + geom_bar(data=df,aes(x = t, y = value, fill = variable),stat = "identity")

# prepare the second layer.
df <- data.frame(t=last(index(v),length_graph),
              value=last(mapply(func,v[,1],v[,2],v[,3],v[,4],v[,5],v[,6],v[,7],v[,8]),length_graph)
)
df <- df[-length(df[,1]),]  # cut off the last entry.
# g <- ggplot(df, aes(x = t, y = value))
g <- g+geom_line(data=df, aes(x = t, y = value))

# plot(g)
png("~/Dropbox/R-script/covid/05tokyo_age.png", width = 1200, height = 800)
plot(g)
dev.off()


2020年6月13日土曜日

EPS 2020JUN13



eps_year_xts["2020::"]
             [,1]
2020-01-01 116.42
2020-04-01 100.84
2020-07-01  93.89
2020-10-01  90.33
2021-01-01 110.84
2021-04-01 125.10
2021-07-01 135.49
2021-10-01 142.90

~$ tac eps.txt 
3/31/2019 2834.40 $37.99 $35.02 18.52 21.09 $153.05 $134.39
6/30/2019 2941.76 $40.14 $34.93 19.04 21.75 $154.54 $135.27
9/30/2019  2976.74 $39.81 $33.99 19.46 22.40 $152.97 $132.90
12/31/2019 3230.78 $39.18 $35.53 20.56 23.16 $157.12 $139.47
3/31/2020 (97.5%) 2584.59 $19.70 $11.98 22.42 26.73 $138.83 $116.43
6/30/2020 $23.08 $19.94 25.56 30.68 $121.77 $101.44
9/30/2020 $31.04 $27.52 27.54 32.77 $113.00 $94.97
12/31/2020 $35.89 $32.62 28.37 33.81 $109.71 $92.06
3/31/2021 $37.02 $33.16 24.50 27.49 $127.03 $113.24
6/30/2021 $38.14 $34.09 21.90 24.43 $142.09 $127.39
9/30/2021 $42.05 $37.98 20.33 22.58 $153.10 $137.85
12/31/2021 $44.59 $40.05 19.24 21.42 $161.80 $145.28

~$   tac eps.txt | awk '{gsub("\\$","",$NF);print "eps_year_xts[\"2019::\"]["NR"] <- "$NF}'

eps_year_xts["2019::"][1] <- 134.39
eps_year_xts["2019::"][2] <- 135.27
eps_year_xts["2019::"][3] <- 132.90
eps_year_xts["2019::"][4] <- 139.47
eps_year_xts["2019::"][5] <- 116.43
eps_year_xts["2019::"][6] <- 101.44
eps_year_xts["2019::"][7] <- 94.97
eps_year_xts["2019::"][8] <- 92.06
eps_year_xts["2019::"][9] <- 113.24
eps_year_xts["2019::"][10] <- 127.39
eps_year_xts["2019::"][11] <- 137.85
eps_year_xts["2019::"][12] <- 145.28

Personal Income and GDP PI



Recalculate based upon the past 20 years except 2020.


summary(lm(last(yearlyReturn(PI),20)[-20] ~ last(yearlyReturn(G),20)[-20]))

Call:
lm(formula = last(yearlyReturn(PI), 20)[-20] ~ last(yearlyReturn(G), 
    20)[-20])

Residuals:
      Min        1Q    Median        3Q       Max 
-0.063473 -0.008838  0.000614  0.013313  0.045431 

Coefficients:
                               Estimate Std. Error t value Pr(>|t|)   
(Intercept)                    0.006467   0.012627   0.512  0.61514   
last(yearlyReturn(G), 20)[-20] 0.878119   0.290504   3.023  0.00767 **
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.02297 on 17 degrees of freedom
Multiple R-squared:  0.3496, Adjusted R-squared:  0.3113 
F-statistic: 9.137 on 1 and 17 DF,  p-value: 0.007673

  • calculate PI according to the formula above on the assumption that 2020 GDPC1 year return is -5% and CPI is 2%.
0.006467 0.878119 *G["2019-10-01"]*0.97
                GDP
2019-10-01 18508.34


  • when 2020 GDP growth is forecasted as -0.05, to calculate PI growth rate do this.


0.006467 + 0.878119 *(expected GDP growth rate)


  • then 2020 PI forecast is as below


(1 + 0.006467 + 0.878119 *-0.05)*PI["2019-10-01"]
                 PI
2019-10-01 18061.11

plot.default(last(yearlyReturn(G),20)[-20],last(yearlyReturn(PI),20)[-20])
abline(lm(last(yearlyReturn(PI),20)[-20] ~ last(yearlyReturn(G),20)[-20]))
abline(v=0)



Recalculate based upon the past 50 years except 2020.


summary(lm(last(yearlyReturn(PI),50)[-50] ~ last(yearlyReturn(G),50)[-50]))

Call:
lm(formula = last(yearlyReturn(PI), 50)[-50] ~ last(yearlyReturn(G), 
    50)[-50])

Residuals:
      Min        1Q    Median        3Q       Max 
-0.065914 -0.005273  0.000005  0.006657  0.043099 

Coefficients:
                               Estimate Std. Error t value Pr(>|t|)    
(Intercept)                    0.008346   0.005670   1.472    0.148    
last(yearlyReturn(G), 50)[-50] 0.890811   0.080517  11.064 1.11e-14 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.01723 on 47 degrees of freedom
Multiple R-squared:  0.7226, Adjusted R-squared:  0.7167 
F-statistic: 122.4 on 1 and 47 DF,  p-value: 1.11e-14

0.008346+0.890811 *G["2019-10-01"]*0.97
                GDP
2019-10-01 18775.85

0.008346+0.890811 *(expected GDP growth rate)
  • when 2020 is -0.05
0.008346+0.890811 *-0.05
[1] -0.03619455

(1+ 0.008346+0.890811 *-0.05)*PI["2019-12-01"]

      PI
2019-12-01 18222.28

2020年6月2日火曜日

.Last function 終了 終了処理



save image files with timestamp file name like "yyyy-mm-dd-hh-mm-ss-00" and .Rdata.

.Last <- function() {
  # save.image(file=paste(getwd(),Sys.time(),sep="/"))
  #
  # save working image in the parent directory. expected the length of the directory name is fixed.
  # adjust value "len" to fit the length of the result "getwd()".
  #
  op <- options(digits.secs = 2)
  start_pos <- 1
  # len <- 23
  # change to decide based upon strings length.
  last_pos <- nchar(as.character(getwd())) -2
  save.image(file=gsub("[ :.]","-",paste(substr(getwd(),start_pos,last_pos),Sys.time(),sep="/")))
   cat("time stamp file done...\n")
  save.image(file=paste(getwd(),".Rdata",sep="/"))
   cat("RData done...\n")
  #
  # OR store working image at the working directory
  #
  # save.image(file=gsub(" ","-",paste(getwd(),Sys.time(),sep="/")))
  #
  cat("bye bye...\n")
}

2020年6月1日月曜日

EPS 2020JUN01



the entries below are wrong.

> eps_year_xts["2020::"]
             [,1]
2020-01-01 102.44
2020-04-01  96.80
2020-07-01  95.16
2020-10-01 115.77
2021-01-01 130.19
2021-04-01 140.44
2021-07-01 147.68
2021-10-01 150.31

~$ tac eps.txt | awk '{gsub("\\$","",$NF);print "eps_year_xts[\"2019::\"]["NR"] <- "$NF}'
eps_year_xts["2019::"][1] <- 134.39
eps_year_xts["2019::"][2] <- 135.27
eps_year_xts["2019::"][3] <- 132.90
eps_year_xts["2019::"][4] <- 139.47
eps_year_xts["2019::"][5] <- 116.42
eps_year_xts["2019::"][6] <- 100.84
eps_year_xts["2019::"][7] <- 93.89
eps_year_xts["2019::"][8] <- 90.33
eps_year_xts["2019::"][9] <- 110.84
eps_year_xts["2019::"][10] <- 125.10
eps_year_xts["2019::"][11] <- 135.49
eps_year_xts["2019::"][12] <- 142.90

2020年5月29日金曜日

Global PMI manufacturing



> pmi_global <- as.xts(c(50.3,50.1,50.4, 47.1,47.3,39.8),seq(as.Date("2019-11-01"),as.Date("2020-04-01"),by='months'))
> pmi_global
           [,1]
2019-11-01 50.3
2019-12-01 50.1
2020-01-01 50.4
2020-02-01 47.1
2020-03-01 47.3
2020-04-01 39.8


> pmi_global <- append(pmi_global,as.xts(c(47.9,50.3,51.8,52.3),seq(as.Date("2020-06-01"),as.Date("2020-09-01"),by='months')))

> pmi_global
           [,1]
2019-11-01 50.3
2019-12-01 50.1
2020-01-01 50.4
2020-02-01 47.1
2020-03-01 47.3
2020-04-01 39.6
2020-05-01 42.4
2020-06-01 47.9
2020-07-01 50.3
2020-08-01 51.8
2020-09-01 52.3