2019年12月17日火曜日

EPS 2019DEC17



> eps_year_xts["2019::"]
             [,1]
2019-01-01 134.39
2019-04-01 135.27
2019-07-01 133.95
2019-10-01 142.58
2020-01-01 146.27
2020-04-01 152.08
2020-07-01 159.78
2020-10-01 165.15

$tac eps.txt 
6/30/2019 2941.76 $40.14 $34.93 19.04 21.75 $154.54 $135.27
9/30/2019 (Prlim.) 2976.74 $39.82 $33.99 20.71 23.84 $152.98 $132.90
12/31/2019 $40.42 $36.51 20.01 22.56 $158.37 $140.45
3/31/2020 $40.60 $36.97 19.68 22.25 $160.98 $142.40
6/30/2020 $43.65 $39.88 19.26 21.50 $164.49 $147.35
9/30/2020 $45.92 $42.64 18.57 20.31 $170.59 $156.01
12/31/2020 $46.35 $42.37 17.95 19.58 $176.52 $161.87

When entries started with "6/30/2019", I have to put NR+1 in order to adjust sequence number.

$tac eps.txt | awk '{gsub("\\$","",$NF);print "eps_year_xts[\"2019::\"]["NR+1"] <- "$NF}'
eps_year_xts["2019::"][2] <- 135.27
eps_year_xts["2019::"][3] <- 132.90
eps_year_xts["2019::"][4] <- 140.45
eps_year_xts["2019::"][5] <- 142.40
eps_year_xts["2019::"][6] <- 147.35
eps_year_xts["2019::"][7] <- 156.01
eps_year_xts["2019::"][8] <- 161.87

2019年12月5日木曜日

xts , matrix and data frame index and dimensions.


# the sample of xts object.

> cli_xts["2019::"]
               oecd      usa    china     ea19
2019-01-01 99.53469 99.71027 98.70411 99.77903
2019-02-01 99.43623 99.54673 98.69403 99.66678
2019-03-01 99.36298 99.40987 98.71894 99.55812
2019-04-01 99.30513 99.28092 98.75597 99.45299
2019-05-01 99.25108 99.15568 98.80080 99.34872
2019-06-01 99.20237 99.03459 98.86234 99.24314
2019-07-01 99.15969 98.91864 98.94160 99.14515
2019-08-01 99.12648 98.81747 99.02364 99.06464
2019-09-01 99.10692 98.75613 99.10512 99.00034

# in order to create the sample matirx, use matrix() to transform xts to marix

> matrix(cli_xts["2019::"],ncol=4)
          [,1]     [,2]     [,3]     [,4]
 [1,] 99.53469 99.71027 98.70411 99.77903
 [2,] 99.43623 99.54673 98.69403 99.66678
 [3,] 99.36298 99.40987 98.71894 99.55812
 [4,] 99.30513 99.28092 98.75597 99.45299
 [5,] 99.25108 99.15568 98.80080 99.34872
 [6,] 99.20237 99.03459 98.86234 99.24314
 [7,] 99.15969 98.91864 98.94160 99.14515
 [8,] 99.12648 98.81747 99.02364 99.06464
 [9,] 99.10692 98.75613 99.10512 99.00034

# matrix require 2 dim. index like [1,] to pick up the designated row.

> matrix(cli_xts["2019::"],ncol=4)[1,]
[1] 99.53469 99.71027 98.70411 99.77903

# while xts requires [1].

> cli_xts["2019::"][1]
               oecd      usa    china     ea19
2019-01-01 99.53469 99.71027 98.70411 99.77903

> cli_xts["2019::"][1,1]
               oecd
2019-01-01 99.53469

> matrix(cli_xts["2019::"],ncol=4)[1,1]
[1] 99.53469

# for the case of data frame.

> data.frame(cli_xts["2019::"],ncol=4)
               oecd      usa    china     ea19 ncol
2019-01-01 99.53469 99.71027 98.70411 99.77903    4
2019-02-01 99.43623 99.54673 98.69403 99.66678    4
2019-03-01 99.36298 99.40987 98.71894 99.55812    4
2019-04-01 99.30513 99.28092 98.75597 99.45299    4
2019-05-01 99.25108 99.15568 98.80080 99.34872    4
2019-06-01 99.20237 99.03459 98.86234 99.24314    4
2019-07-01 99.15969 98.91864 98.94160 99.14515    4
2019-08-01 99.12648 98.81747 99.02364 99.06464    4
2019-09-01 99.10692 98.75613 99.10512 99.00034    4

# a index works in the same way as matrix.

> data.frame(cli_xts["2019::"],ncol=4)[1,]
               oecd      usa    china     ea19 ncol
2019-01-01 99.53469 99.71027 98.70411 99.77903    4


> data.frame(cli_xts["2019::"],ncol=4)[1,1]
[1] 99.53469

# a single dimention index also works fine for data frame
# but it takes column for the parameter

> data.frame(cli_xts["2019::"],ncol=4)[1]
               oecd
2019-01-01 99.53469
2019-02-01 99.43623
2019-03-01 99.36298
2019-04-01 99.30513
2019-05-01 99.25108
2019-06-01 99.20237
2019-07-01 99.15969
2019-08-01 99.12648
2019-09-01 99.10692

# on the other hand, single dim. index taken for the row.

> cli_xts["2019::"][1]
               oecd      usa    china     ea19
2019-01-01 99.53469 99.71027 98.70411 99.77903

2019年11月14日木曜日

序数の算出 how to find sequence number from data frame


wpm <- rbind(merge(wp,rep(1,length(index(wp))),suffixes =c('','p_or_m')),merge(wm,rep(-1,length(index(wm)))))
wpm
           result period_length          rate  open_p close_p p_or_m
 9 1970 0.9157072             9 -0.0097365802   92.06   84.30     -1
 1 1973 1.3763938            28  0.0114748680   84.30  116.03      1
 1 1975 0.6634491            24 -0.0169506551  116.03   76.98     -1
 5 1976 1.3013769            16  0.0166002074   76.98  100.18      1
 7 1977 0.9867239            14 -0.0009541901  100.18   98.85     -1
                    <skip>
 4 2016 1.0031085            16  0.0001940004 2058.90 2065.30     -1
12 2017 1.2933674            20  0.0129455360 2067.17 2673.61      1
 8 2019 1.0904450            19  0.0045675511 2683.73 2926.46     -1

# return 1.0 以上のマイナスフェーズを探す。
# return 1.0 以上のマイナスフェーズを探して、その日付を取得する。

 index(wpm[,1])[wpm[,6] == -1 & wpm[,1] > 1]

 [1] " 6 1980" "11 1984" " 7 1986" "10 1989" " 1 1991" "12 1991" "11 1992" " 6 1995" " 1 1996" "10 1998" " 5 2005" " 6 2006"
[13] " 7 2010" " 9 2012" " 8 2014" " 4 2016" " 8 2019"

# return 1.0 以上のマイナスフェーズを探して、その日付ではなく
# その序数を取得する。51が最終エントリを指す序数であることに注意。

 index(as.vector(wpm[,1]))[wpm[,6] == -1 & wpm[,1] > 1]

 [1]  7 11 13 17 19 21 23 25 27 29 35 37 41 45 47 49 51

# 取得した序数に+1し、その次のエントリの序数を得る。

  index(as.vector(wpm[,1]))[wpm[,6] == -1 & wpm[,1] > 1]+1

 [1]  8 12 14 18 20 22 24 26 28 30 36 38 42 46 48 50 52

#シーケンスの要素数を調べる。51が最終エントリなので、序数=51+1に該当する要素は
# 元のデータフレームにはないことに注意。

  length(index(as.vector(wpm[,1]))[wpm[,6] == -1 & wpm[,1] > 1]+1)
[1] 17

#[-17]をつけて該当する序数を削除の上、データを再構成する

wpm[(index(as.vector(wpm[,1]))[wpm[,6] == -1 & wpm[,1] > 1]+1)[-17],]

           result period_length         rate  open_p close_p p_or_m
12 1980 1.1883753             6  0.029182211  114.24  135.76      1
 2 1985 1.1075926             3  0.034649717  163.58  181.18      1
 8 1987 1.3967474            13  0.026036742  236.12  329.80      1
 2 1990 0.9751147             4 -0.006280247  340.36  331.89      1
 8 1991 1.1498066             7  0.020142134  343.91  395.43      1
 4 1992 0.9950124             4 -0.001249244  417.03  414.95      1
 9 1994 1.0727019            22  0.003195123  431.35  462.71      1
10 1995 1.0674621             4  0.016454915  544.75  581.50      1
 9 1997 1.4893871            20  0.020117927  636.02  947.28      1
 3 2000 1.3639946            17  0.018427587 1098.67 1498.58      1
 5 2006 1.0659588            12  0.005337085 1191.50 1270.09      1
 6 2007 1.1836842            12  0.014151848 1270.06 1503.35      1
 2 2011 1.1983603             7  0.026187620 1107.53 1327.22      1
 3 2014 1.2994239            18  0.014657552 1440.90 1872.34      1
12 2014 1.0273593             4  0.006770750 2004.07 2058.90      1
12 2017 1.2933674            20  0.012945536 2067.17 2673.61      1

Download composite leading indicator full data and load them into the R environment CLI

# go to OECD, and download 'Full Indicator data(.csv)".
# change the filename to CLI3.csv upon download.
# go to the directory where the downloaded csv exists.
# run the script in git repo. for example as the below.

$~/R/R2/index/cli_download.sed

# there will be

$ls -l *.csv
-rw-r--r--@ 1 honomoto  staff  1283086 11 14 09:50 CLI3.csv
-rw-r--r--  1 honomoto  staff     7105 11 14 09:50 chn.csv
-rw-r--r--  1 honomoto  staff    12848 11 14 09:50 ea19.csv
-rw-r--r--  1 honomoto  staff    13923 11 14 09:50 oecd.csv
-rw-r--r--  1 honomoto  staff    15465 11 14 09:50 usa.csv

# run cli_download.


cli_xts <- merge(as.xts(read.zoo(read.csv("~/Downloads/oecd.csv"))),
as.xts(read.zoo(read.csv("~/Downloads/usa.csv"))),
as.xts(read.zoo(read.csv("~/Downloads/chn.csv"))),
as.xts(read.zoo(read.csv("~/Downloads/ea19.csv"))),
suffixes = c("oecd","usa","china","ea19"))

# draw graphs

par(mfrow=c(4,1))
plot(diff(cli_xts$oecd)["2011::"],type='p',pwd=1,pch='+')
plot(diff(cli_xts$usa)["2011::"],type='p',pwd=1,pch='+')
plot(diff(cli_xts$ea19)["2011::"],type='p',pwd=1,pch='+')
plot(diff(cli_xts$china)["2011::"],type='p',pwd=1,pch='+')

par(mfrow=c(1,1))

2019年11月13日水曜日

Download Composite Leading Indicator from OECD site. Early Bird Version


# go to OECD site and download CSV from Export menu.
# assume download file name is  "MEI_CLI.csv".
# pick up amplitute justfied OECD total data from CSV

$sed -n '/OECD\ -\ Total/p' MEI_CLI_13112019035737453.csv |grep LOLITOAA | grep Ampli | awk -F, 'BEGIN{ORS = ""}{print $(NF-2)","}'

the improved version is below to integrate all query string into one criteria.

sed -n '/LOLITOAA.*Ampli.*OECD\ -\ Total/p' MEI_CLI.csv | awk -F, 'BEGIN{ORS = "";print "c("}{print $(NF-2)","}END{print $(NF-2)")\n"}'


OR do as below.



awk -F, 'BEGIN{ORS = "";print "w <- c("}/LOLITOAA.*Ampli.*OECD\ -\ Total/{print $(NF-2)","}END{print $(NF-2)")\n length(w)\n"}' MEI_CLI.csv 

output sample is like below. now output is in the from of R statement. Construct vector in c() and substitute to w().


100.7565,100.7652,100.7551,100.7274,100.6736,100.5974,100.512,100.4191,100.3163,100.2006,100.0703,99.93036,99.79047,99.6553,99.53469,99.43623,99.36298,99.30513,99.25108,99.20237,99.15969,99.12648,99.10692

w <- c(100.7565,100.7652,100.7551,100.7274,100.6736,100.5974,100.512,100.4191,100.3163,100.2006,100.0703,99.93036,99.79047,99.6553,99.53469,99.43623,99.36298,99.30513,99.25108,99.20237,99.15969,99.12648,99.10692)

length(w)
[1] 23

# The last month will be found in this way.

sed -n '/LOLITOAA.*Ampli.*OECD\ -\ Total/p' MEI_CLI.csv  | awk -F, 'END{print $7}'
"2019-09"

# as the last entry is "2019-09-01" do as below

w <- as.xts(w,last(seq(as.Date("2010-01-01"),as.Date("2019-09-01"),by='months'),length(w))

2019年11月12日火曜日

CLI1ヶ月変化値、X軸はCLI5ヶ月変化値、月間騰落率、COV





# 0) Dropbox/R-script/monthreturn_delta5month_delta1monthcolor.r 参照。
# 1)Y軸はCLI1ヶ月変化値、X軸はCLI5ヶ月変化値、同じく月間騰落率を色と形で表す。
# 2)原点からの角度をX、月間騰落率をYに取ったグラフ。
# 3)同じく、2)原点からの角度をX、COVをYに取ったグラフ。
#

w <- (to.monthly(SP5)[,4]/to.monthly(SP5)[,1])["1970::2018"]
c <- (apply.monthly(SP5[,4],sd) / apply.monthly(SP5[,4],mean) )["1970::2018"]
w <- w-1
# w <- (apply.monthly(SP5[,4],sd)/apply.monthly(SP5[,4],mean))["1970::2018"]
d <- na.omit(diff(cli_xts$oecd,5))["1970::2018"]
func <- function(x){
  if(x > 0.1){return("a")}
  if(x > 0.025){return("b")}
  if(x > 0){return("c")}
  if(x > -0.025){return("d")}
  if(x > -0.05){return("e")}
  if(x < -0.05){return("f")}
}
# df <- data.frame(monthlyreturn=as.vector(w),delta=as.vector(d),sign=as.vector(apply(diff(cli_xts$oecd)["1970::2018"],1,func)))

df <- data.frame(monthlyreturn=as.vector(apply(w,1,func)),
          delta=as.vector(d),sign=as.vector(diff(cli_xts$oecd)["1970::2018"]))
  # as.vector(w),delta=as.vector(d),sign=as.vector(apply(diff(cli_xts$oecd)["1970::2018"],1,func)))
p <- ggplot(df, aes(x=delta,y=sign))
p <- p + theme(panel.background = element_rect(fill = "black",
                                               colour = "lightblue"),
               legend.key = element_rect(fill='black',colour='white'))
               # legend.box.background=element_rect(fill='black',colour='lightblue'))

p <- p + geom_point(alpha=1,aes(color=monthlyreturn,shape=monthlyreturn))
p <- p + geom_vline(xintercept =as.vector(last(diff(cli_xts$oecd,5))),size=0.5,linetype=2,colour="white",alpha=0.5)
p <- p + geom_hline(yintercept =as.vector(last(diff(cli_xts$oecd))),size=0.5,linetype=2,colour="white",alpha=0.5)
p <- p + theme(axis.title.x=element_blank(),axis.title.y=element_blank())
# p <- p + geom_smooth(method = "auto")
# p <- p + scale_shape_manual(values=c(0,1,2,10,11,12))
# p <- p + scale_color_brewer(label=c("more than 0.1","more than 0.025","more than ZERO","more then -0.025","more than -0.1","less than -0.1"))  #x-axis label
p <- p + scale_color_brewer(palette="Spectral", label=c("more than 0.1","more than 0.025","more than ZERO","more then -0.025","more than -0.05","less than -0.05"))  #x-axis label
p <- p + scale_shape_manual(values=c(0,1,2,10,11,12),label=c("more than 0.1","more than 0.025","more than ZERO","more then -0.025","more than -0.05","less than -0.05"))
p <- p + stat_smooth(aes(x=delta,y=sign),method="loess",color='white',size=0.3,se=FALSE)
# p <- p + guides(shape = FALSE)
plot(p)




wdf <- data.frame(dist=sqrt(df$delta**2+df$sign**2),
                  angle=atan2(df$delta,df$sign),
                  return=w,
                  mon=seq(1,length(wdf[,1]),1),
                  cov=c)
colnames(wdf)[3] <- "return"
# wdf <- cbind(wdf,seq(1,length(wdf[,1]),1))
colnames(wdf) [4] <- "mon"
# wdf <- cbind(wdf,c)
# wdf
colnames(wdf)[5] <- "cov"


df <- data.frame(r<-wdf[,3],a<-wdf[,2],mon<-wdf[,4],cov <- wdf[,5])
colnames(df)[1] <- 'r'
colnames(df)[2] <- 'a'
colnames(df)[3] <- 'mon'
colnames(df)[4] <- 'c'

# df <- last(df,284)
p <- ggplot(df,aes(x=a,y=r))
p <- p + geom_point(alpha=1,aes(color=mon))
p <- p + geom_vline(xintercept =as.vector(atan2(last(diff(cli_xts$oecd,5),1),last(diff(cli_xts$oecd),1))),size=0.5,linetype=2,colour="red",alpha=0.5)
p <- p + geom_vline(xintercept =as.vector(seq(-0.9,0.9,0.1))*pi,size=0.5,linetype=2,colour="white",alpha=0.5)
# p <- p + scale_color_brewer(palette="Spectral")
p <- p + stat_smooth(aes(x=a,y=r),method="loess",color='white',size=0.3,se=FALSE)
p <- p + scale_color_gradient(low = "red", high = "green")
p <- p + theme(panel.background = element_rect(fill = "black",
                                              colour = "lightblue"))
             # panel.grid = element_blank())
plot(p)

p <- ggplot(df,aes(x=a,y=c))
p <- p + geom_point(alpha=1,aes(color=mon))
p <- p + geom_vline(xintercept =as.vector(atan2(last(diff(cli_xts$oecd,5),1),last(diff(cli_xts$oecd),1))),size=0.5,linetype=2,colour="red",alpha=0.5)
p <- p + geom_vline(xintercept =as.vector(seq(-0.9,0.9,0.1))*pi,size=0.5,linetype=2,colour="white",alpha=0.5)
# p <- p + scale_color_brewer(palette="Spectral")
p <- p + stat_smooth(aes(x=a,y=c),method="loess",color='white',size=0.3,se=FALSE)
p <- p + scale_color_gradient(low = "red", high = "green")
p <- p + theme(panel.background = element_rect(fill = "black",
                                              colour = "lightblue"))
             # panel.grid = element_blank())
plot(p)

2019年11月11日月曜日

EPS 2019NOV11


before

> eps_year_xts["2019::"]
             [,1]
2019-01-01 134.39
2019-04-01 135.27
2019-07-01 135.99
2019-10-01 145.71
2020-01-01 148.71
2020-04-01 154.31
2020-07-01 159.83
2020-10-01 165.38

# Edit delete 2 lines in the middle


12/31/2020 $46.68 $42.96 17.40 18.68 $177.27 $165.15
9/30/2020 $45.92 $42.75 18.01 19.31 $171.35 $159.78
6/30/2020 $43.85 $40.74 18.62 20.29 $165.72 $152.08
3/31/2020 $40.82 $38.70 19.04 21.09 $162.01 $146.27
12/31/2019 $40.76 $37.59 19.38 21.64 $159.18 $142.58
9/30/2019 (89.2%) 2976.74 $40.29 $35.04 20.11 23.03 $153.45 $133.95
19.40 22.22 (P/E on Sep,'19 price)
ACTUALS
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


# result is like this.

honomoto@~$cat eps.txt
12/31/2020 $46.68 $42.96 17.40 18.68 $177.27$165.15
9/30/2020 $45.92 $42.75 18.01 19.31 $171.35$159.78
6/30/2020 $43.85 $40.74 18.62 20.29 $165.72$152.08
3/31/2020 $40.82 $38.70 19.04 21.09 $162.01$146.27
12/31/2019 $40.76 $37.59 19.38 21.64 $159.18$142.58
9/30/2019 (89.2%) 2976.74 $40.29 $35.04 20.11 23.03 $153.45 $133.95
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




honomoto@~$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] <- 133.95
eps_year_xts["2019::"][4] <- 142.58
eps_year_xts["2019::"][5] <- 146.27
eps_year_xts["2019::"][6] <- 152.08
eps_year_xts["2019::"][7] <- 159.78
eps_year_xts["2019::"][8] <- 165.15


2019年11月6日水曜日

確率密度の計算 dnorm() qnorm()


オブジェクト wdf はこちらから。

以下のようなデータの集合があったとする。また、そのヒストグラムを示す。

> wdf[,3][wdf[,2] > -0.5*pi & wdf[,2] < 1]

[17]  0.0711926070 -0.0854386059  0.0094247436  0.0030609850  0.0432560109  0.0165413246  0.0213577174  0.0673432082
[25]  0.0222028328  0.0003197991 -0.0197827495  0.0095648204 -0.0217844120  0.0101077662  0.0070690542  0.0104836274
[33]  0.0317760404 -0.0003202381  0.0400974724  0.0069337016  0.0079165671  0.0134314222  0.0591259963  0.0563752866
[41]  0.0410094415  0.0180991114  0.0751759043  0.0075738342  0.0810441069  0.0508987097  0.0113222145 -0.0001427142
[49]  0.0359682872 -0.0112221049  0.0936167706  0.0998248066  0.0531659206  0.0085731462  0.0425385266 -0.0199458599
[57]  0.0028183707  0.0069545273 -0.0158576835  0.0236426721  0.0244524103  0.0144061834  0.0023497169  0.0353730253


その時平均と分散はそれぞれmean()とsd()によって算出される。

> mean(wdf[,3][wdf[,2] > -0.5*pi & wdf[,2] < 1])
[1] 0.02359295
> sd(wdf[,3][wdf[,2] > -0.5*pi & wdf[,2] < 1])
[1] 0.03504609

念のため、シャピロウイルク検定を実施する。このテストの帰無仮説は「正規分布は成立している」だから注意すること。結果は9.5%の確率で正規分布が成立することになる。あまり良い値ではない。

> shapiro.test(wdf[,3][wdf[,2] > -0.5*pi & wdf[,2] < 1])

Shapiro-Wilk normality test

data:  wdf[, 3][wdf[, 2] > -0.5 * pi & wdf[, 2] < 1]
W = 0.96799, p-value = 0.095



平均と分散が判明すればそれに応じて確率を計算することが可能になる。

 下方累積確率 = pnorm(標本値、平均、分散)

> pnorm(0,mean(wdf[,3][wdf[,2] > -0.5*pi & wdf[,2] < 1]), sd(wdf[,3][wdf[,2] > -0.5*pi & wdf[,2] < 1]))
[1] 0.2504107

正規分布に従えば、平均と分散が上に示された集団で標本の値が0以下となる確率は25.04%となる。

> pnorm(0.02359295,mean(wdf[,3][wdf[,2] > -0.5*pi & wdf[,2] < 1]), sd(wdf[,3][wdf[,2] > -0.5*pi & wdf[,2] < 1]))
[1] 0.5

また、値が0.02359295以下となる確率は50%である。また逆に確率からありうべき標本の値を計算することも可能になる。確率50%を期待できるのは、無限小から0.02359295までの累積となる。

想定標本値= qnorm(下方累積確率、平均、分散)

> qnorm(0.5,mean(wdf[,3][wdf[,2] > -0.5*pi & wdf[,2] < 1]), sd(wdf[,3][wdf[,2] > -0.5*pi & wdf[,2] < 1]))
[1] 0.02359295

2019年11月5日火曜日

rolling COV , grit.ticks.on='months', sapply



#
# calculate OCV of each last 30 business days with given parameters.
#
func <- function(date,xts=SP5[,4],days=30){
  hensa <- sd(last(xts[paste("::",date,sep="")],days))
  heikin <- mean(last(xts[paste("::",date,sep="")],days))
  return(hensa/heikin)
}
#
# construct XTS object
# use sapply() instead of for loop to calculate rolling COV
#
as.xts(sapply(last(index(SP5),365),func),last(index(SP5),365))
# plot!!
plot(as.xts(sapply(last(index(SP5),730),func),last(index(SP5),730)),main="moving COV",grid.ticks.on="months")
length(index(SP5["2012-01-01::"]))
#
# 1)automatically the number of business days from the specific date (2012/12/1).
# 2)calculate rolling COV from the date above.
# 3)plot graph with monthly grid.
#
plot(as.xts(sapply(last(index(SP5),length(index(SP5["2012-01-01::"]))),func),last(index(SP5),length(index(SP5["2012-01-01::"])))),main="moving COV",grid.ticks.on="months")

# 評価比較用 for を使用した場合
system.time(for( i in as.character(last(index(SP5),365))){cat(i); print(func(i))})

system.time(print(as.xts(sapply(last(index(SP5),365),func),last(index(SP5),365))))

2019年10月30日水曜日

calculate S&P 500 n month return.

w <- c()
for(i in seq(1,13,1)){
 m <- matrix(as.vector(to.monthly(SP5)[,4])[i:(i+818)],nrow=13,ncol=63)
  w <- append(
    w,m[13,] / m[1,]

  )
}


gap <- 61 # when month gap is 61.
a<-Sys.time()
w <- c()
j <- (floor(length(index(to.monthly(SP5)))/gap)-1)*gap-1
for(i in seq(1,gap,1)){
#  m <- matrix(as.vector(to.monthly(SP5)[,4])[i:(i+j)],nrow=gap,ncol=(floor(length(index(to.monthly(SP5)))/gap)-1))
m <- matrix(as.vector(to.monthly(SP5)[,4])[i:(i+j)],nrow=gap,ncol=((j+1)/gap))
  w <- append(
    w,m[gap,] / m[1,]
  )
}

s <- as.vector(matrix(t(matrix(w,ncol=gap)),ncol=1))

length(s[s<1])
mean(s[s<1])
mean(s[s>1])
length(s[s>1])
b<-Sys.time()
b -a

# テスト用サンプル gap = 61のとき有効。
# > s[1:60]
#  [1] 2.148387 2.134727 2.115674 2.113586 2.018637 2.319389 2.439462 2.344191 2.245244 2.167947
# [11] 2.332650 2.226138 2.023084 2.079817 2.256983 2.156933 2.100372 2.240935 2.204911 2.040807
# [21] 1.949699 1.986922 1.970280 1.963399 1.852527 1.859845 1.810012 1.961407 1.987846 1.897837
# [31] 1.886220 1.806632 1.728606 1.674551 1.625877 1.505081 1.580743 1.576834 1.664689 1.764419
# [41] 1.796658 1.874068 1.906667 2.047599 2.143897 2.091687 2.119548 2.225312 2.126150 2.118929
# [51] 2.057906 2.037863 2.010277 2.001712 1.959521 1.997989 1.760446 1.815657 1.702103 1.664536
# > as.vector(to.monthly(SP5)[,4])[61:120]/as.vector(to.monthly(SP5)[,4])[1:60]
#  [1] 2.148387 2.134727 2.115674 2.113586 2.018637 2.319389 2.439462 2.344191 2.245244 2.167947
# [11] 2.332650 2.226138 2.023084 2.079817 2.256983 2.156933 2.100372 2.240935 2.204911 2.040807
# [21] 1.949699 1.986922 1.970280 1.963399 1.852527 1.859845 1.810012 1.961407 1.987846 1.897837
# [31] 1.886220 1.806632 1.728606 1.674551 1.625877 1.505081 1.580743 1.576834 1.664689 1.764419
# [41] 1.796658 1.874068 1.906667 2.047599 2.143897 2.091687 2.119548 2.225312 2.126150 2.118929
# [51] 2.057906 2.037863 2.010277 2.001712 1.959521 1.997989 1.760446 1.815657 1.702103 1.664536

# > s[101:160]
#  [1] 1.605806 1.533378 1.464929 1.518325 1.432281 1.441847 1.395389 1.358812 1.389360 1.404079
# [11] 1.424603 1.379753 1.369632 1.397127 1.374649 1.372987 1.479958 1.475313 1.448524 1.415094
# [21] 1.574537 1.557912 1.556921 1.638955 1.583736 1.477864 1.535759 1.530372 1.680867 1.731036
# [31] 1.649442 1.590604 1.503399 1.437894 1.371503 1.394273 1.294020 1.310953 1.252247 1.132658
# [41] 1.147310 1.168755 1.128015 1.122711 1.258135 1.240423 1.296909 1.440987 1.493879 1.655525
# [51] 1.627168 1.583897 1.718678 1.650743 1.509798 1.528843 1.393354 1.389952 1.354965 1.396275
# > as.vector(to.monthly(SP5)[,4])[161:220]/as.vector(to.monthly(SP5)[,4])[101:160]
#  [1] 1.605806 1.533378 1.464929 1.518325 1.432281 1.441847 1.395389 1.358812 1.389360 1.404079
# [11] 1.424603 1.379753 1.369632 1.397127 1.374649 1.372987 1.479958 1.475313 1.448524 1.415094
# [21] 1.574537 1.557912 1.556921 1.638955 1.583736 1.477864 1.535759 1.530372 1.680867 1.731036
# [31] 1.649442 1.590604 1.503399 1.437894 1.371503 1.394273 1.294020 1.310953 1.252247 1.132658
# [41] 1.147310 1.168755 1.128015 1.122711 1.258135 1.240423 1.296909 1.440987 1.493879 1.655525
# [51] 1.627168 1.583897 1.718678 1.650743 1.509798 1.528843 1.393354 1.389952 1.354965 1.396275

when gap = 61,

> length(s[s<1])
[1] 137
> mean(s[s<1])
[1] 0.8869793
> mean(s[s>1])
[1] 1.629065
> length(s[s>1])
[1] 595

when gap = 14,

> length(s[s<1])
[1] 216
> mean(s[s<1])
[1] 0.8863247
> mean(s[s>1])
[1] 1.160806
> length(s[s>1])
[1] 602
> b<-Sys.time()

2019年10月28日月曜日

CLI delta plus period graph.


# draw wp graph starts.


wpx <- wp

wp <- wp[wp[,2] > 10]

par(mfrow=c(2,3))
plot(to.monthly(SP5[,4])[,4][last(paste(  substr(mondate(index(wp)) - as.vector(wp[,2]-1),1,7), substr(mondate(index(wp)),1,7),sep="::"),1)[1]],main="")
plot(to.monthly(SP5[,4])[,4][last(paste(  substr(mondate(index(wp)) - as.vector(wp[,2]-1),1,7), substr(mondate(index(wp)),1,7),sep="::"),2)[1]],main="")
plot(to.monthly(SP5[,4])[,4][last(paste(  substr(mondate(index(wp)) - as.vector(wp[,2]-1),1,7), substr(mondate(index(wp)),1,7),sep="::"),3)[1]],main="")
plot(to.monthly(SP5[,4])[,4][last(paste(  substr(mondate(index(wp)) - as.vector(wp[,2]-1),1,7), substr(mondate(index(wp)),1,7),sep="::"),4)[1]],main="")
plot(to.monthly(SP5[,4])[,4][last(paste(  substr(mondate(index(wp)) - as.vector(wp[,2]-1),1,7), substr(mondate(index(wp)),1,7),sep="::"),5)[1]],main="")
plot(to.monthly(SP5[,4])[,4][last(paste(  substr(mondate(index(wp)) - as.vector(wp[,2]-1),1,7), substr(mondate(index(wp)),1,7),sep="::"),6)[1]],main="")




plot(to.monthly(SP5[,4])[,4][last(paste(  substr(mondate(index(wp)) - as.vector(wp[,2]-1),1,7), substr(mondate(index(wp)),1,7),sep="::"),7)[1]],main="")
plot(to.monthly(SP5[,4])[,4][last(paste(  substr(mondate(index(wp)) - as.vector(wp[,2]-1),1,7), substr(mondate(index(wp)),1,7),sep="::"),8)[1]],main="")
plot(to.monthly(SP5[,4])[,4][last(paste(  substr(mondate(index(wp)) - as.vector(wp[,2]-1),1,7), substr(mondate(index(wp)),1,7),sep="::"),9)[1]],main="")
plot(to.monthly(SP5[,4])[,4][last(paste(  substr(mondate(index(wp)) - as.vector(wp[,2]-1),1,7), substr(mondate(index(wp)),1,7),sep="::"),10)[1]],main="")
plot(to.monthly(SP5[,4])[,4][last(paste(  substr(mondate(index(wp)) - as.vector(wp[,2]-1),1,7), substr(mondate(index(wp)),1,7),sep="::"),11)[1]],main="")
plot(to.monthly(SP5[,4])[,4][last(paste(  substr(mondate(index(wp)) - as.vector(wp[,2]-1),1,7), substr(mondate(index(wp)),1,7),sep="::"),12)[1]],main="")




plot(to.monthly(SP5[,4])[,4][last(paste(  substr(mondate(index(wp)) - as.vector(wp[,2]-1),1,7), substr(mondate(index(wp)),1,7),sep="::"),13)[1]],main="")
plot(to.monthly(SP5[,4])[,4][last(paste(  substr(mondate(index(wp)) - as.vector(wp[,2]-1),1,7), substr(mondate(index(wp)),1,7),sep="::"),14)[1]],main="")
# plot(to.monthly(SP5[,4])[,4][last(paste(  substr(mondate(index(wp)) - as.vector(wp[,2]-1),1,7), substr(mondate(index(wp)),1,7),sep="::"),15)[1]],main="")
# plot(to.monthly(SP5[,4])[,4][last(paste(  substr(mondate(index(wp)) - as.vector(wp[,2]-1),1,7), substr(mondate(index(wp)),1,7),sep="::"),16)[1]],main="")
# plot(to.monthly(SP5[,4])[,4][last(paste(  substr(mondate(index(wp)) - as.vector(wp[,2]-1),1,7), substr(mondate(index(wp)),1,7),sep="::"),17)[1]],main="")
# plot(to.monthly(SP5[,4])[,4][last(paste(  substr(mondate(index(wp)) - as.vector(wp[,2]-1),1,7), substr(mondate(index(wp)),1,7),sep="::"),18)[1]],main="")


par(mfrow=c(1,1))
wp <- wpx
remove(wpx)


# draw wp graph ends.






2019年10月24日木曜日

総力戦

異次元緩和はともかく、ここ10年ほど世界中で資産の価格は上がるが労働者の賃金は全然上がらない。労働と資本の間のパワーバランスが崩れて、資本が決定的な優位を占める状態となっている。(続

総力戦体制が福祉国家を要求したが、一旦は実現された福祉国家自体が持続可能な国家モデルとしての妥当性を欠いていることが明らかになると同時に、総力戦は時代遅れで不可能になった。今の時代に誰が総力戦を戦うことを真剣に考えているだろうか?(続

同時に、技術進歩は経済を資本集約的なものに変化させた。総力戦に必須だった忠誠心も持ち熟練した技能を持つ大量の労働者はもはや必要とされていない。練度の低い、規格化した労働を、その時々の要求に応じてこなすだけの大衆の群れが残った。(続

今でも必要とされる熟練労働者は高度な教育を受けた知識労働者だ。しかし、彼らは本質的には一種の資本財である。知識集約型産業は形を変えた資本集約型産業に過ぎない。現代の知識労働者は付加価値生産のために莫大な設備投資を必要とする。コンピュータを持たないプログラマーになんの価値があるか?(続

知識と設備、両者のどちらがより本質的かという問いは極めてトリビアルなものと言わざるを得ない。かくして、経済の焦点が労働から資本に移動した以上、資産の価格は上がるが労働者の賃金は全然上がらないのは自然なことである。(続

これらはいずれも本質的な変化なので、所得税の累進割合を変えたり、資産課税を導入したり小手先の政策をいくらいじっても事態の本質は変わらない。今更総力戦を戦いたい人間はいない。しかし、総力戦なき福祉国家の漂流を止めることも誰にもできない。(終

2019年10月8日火曜日

MOV and VIX on each months.



# 各月のCOVを計算して棒グラフにする。
# COVのデータを遡って、12で割り切れる最長のデータ長を計算する。
# データ長に基づいて、配列を作成し、それを数=12の行列に代入する。(行数=12でなくてはならない。注意!)
# apply(行列,1,mean)を実行して、月ごとの平均を計算する。
# 最終COVの月を計算して、データの並びを変える。最終月が10月なら、3,4,5,6,7,8,9,10,11,12,1,2となるように並べる。



df <- data.frame(cov=apply(matrix(last(COV,floor(length(COV)/12)*12),nrow=12),1,mean)[c(last(seq(1,12,1),month(last(index(COV)))),seq(1,(12 - month(last(index(COV)))),1))],
month=factor(seq(1,12,1)))
p <- ggplot(df,aes(x=month,fill=month))
p <- p + geom_bar(aes(y=cov),stat='identity')
p <- p + scale_x_discrete(label=c("jan","feb","mar","apr","may","jun","jul","aug","sep","oct","nov","dec"))  #x-axis label
p <- p + theme(legend.position = 'none')  # erase legend
plot(p)





df <- data.frame(cov=apply(matrix(last(VIX[,4],floor(length(VIX[,4])/12)*12),nrowncol=12),12,mean)[c(last(seq(1,12,1),month(last(index(VIX[,4])))),seq(1,(12 - month(last(index(VIX[,4])))),1))],
month=factor(seq(1,12,1)))
p <- ggplot(df,aes(x=month,fill=month))
p <- p + geom_bar(aes(y=cov),stat='identity')
p <- p + scale_x_discrete(label=c("jan","feb","mar","apr","may","jun","jul","aug","sep","oct","nov","dec"))  #x-axis label
p <- p + theme(legend.position = 'none')  # erase legend
plot(p)



2019年10月7日月曜日

EPS 2019OCT07

$cat eps.txt 

12/31/2020 $47.81 $44.23 16.53 18.00 $180.06 $165.38
9/30/2020 $46.25 $42.60 17.06 18.62 $174.46 $159.83
6/30/2020 $44.53 $40.53 17.61 19.29 $169.02 $154.31
3/31/2020 $41.47 $38.02 18.08 20.02 $164.63 $148.71
12/31/2019 $42.21 $38.69 18.47 20.43 $161.15 $145.71
9/30/2019 2976.74 $40.81 $37.08 19.33 21.89 $153.97 $135.99
6/30/2019 (97.7%) 2941.76 $40.14 $34.93 19.26 22.01 $154.54 $135.27

3/31/2019 2834.40 $37.99 $35.02 18.76 21.30 $153.05 $134.39

# use tac command to invert eps.txt.



$ tac eps.txt 

3/31/2019 2834.40 $37.99 $35.02 18.76 21.30 $153.05 $134.39
6/30/2019 (97.7%) 2941.76 $40.14 $34.93 19.26 22.01 $154.54 $135.27
9/30/2019 2976.74 $40.81 $37.08 19.33 21.89 $153.97 $135.99
12/31/2019 $42.21 $38.69 18.47 20.43 $161.15 $145.71
3/31/2020 $41.47 $38.02 18.08 20.02 $164.63 $148.71
6/30/2020 $44.53 $40.53 17.61 19.29 $169.02 $154.31
9/30/2020 $46.25 $42.60 17.06 18.62 $174.46 $159.83
12/31/2020 $47.81 $44.23 16.53 18.00 $180.06 $165.38

$ 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] <- 135.99
eps_year_xts["2019::"][4] <- 145.71
eps_year_xts["2019::"][5] <- 148.71
eps_year_xts["2019::"][6] <- 154.31
eps_year_xts["2019::"][7] <- 159.83
eps_year_xts["2019::"][8] <- 165.38


> eps_year_xts["2019::"]
             [,1]
2019-01-01 134.39
2019-04-01 135.44
2019-07-01 136.41
2019-10-01 146.35
2020-01-01 149.59
2020-04-01 155.21
2020-07-01 160.75
2020-10-01 166.45
> eps_year_xts["2019::"][1] <- 134.39
> eps_year_xts["2019::"][2] <- 135.27
> eps_year_xts["2019::"][3] <- 135.99
> eps_year_xts["2019::"][4] <- 145.71
> eps_year_xts["2019::"][5] <- 148.71
> eps_year_xts["2019::"][6] <- 154.31
> eps_year_xts["2019::"][7] <- 159.83
> eps_year_xts["2019::"][8] <- 165.38
> eps_year_xts["2019::"]
             [,1]
2019-01-01 134.39
2019-04-01 135.27
2019-07-01 135.99
2019-10-01 145.71
2020-01-01 148.71
2020-04-01 154.31
2020-07-01 159.83
2020-10-01 165.38

2019年9月17日火曜日

SPX from 2000 + residual with eps model price + CLI 1 month delta.


SPX from 2000 + residual with eps model price + CLI 1 month delta.

kikan <- "2000-01-01::"
func <- function(x){
  if(is.na(x)){return(NA)}
  if(x > 0.1){return("a")}
  if(x > 0){return("b")}
  if(x > -0.1){return("c")}
  if(x < -0.1){return("d")}
}

delta <- append(as.vector(diff(cli_xts$oecd)[kikan]),rep(NA,length(index(tmp.predict)) - length(diff(cli_xts$oecd)[kikan])))
#  [kikan]
df <- data.frame(i=as.vector(tmp.predict[kikan][,4]),
g=as.vector(tmp.predict[kikan][,6]),
e=as.vector(tmp.predict[kikan][,7]),
r=as.vector((tmp.predict[kikan][,4]/tmp.predict[kikan][,7])-1)*5000,
# c=as.vector(apply(diff(cli_xts$oecd)[kikan],1,func)),
c=as.vector(apply(matrix(delta,ncol=1),1,func)),
t=as.Date(index(tmp.predict[kikan])))

colnames(df)[1] <- 'i'
colnames(df)[2] <- 'g'
colnames(df)[3] <- 'e'
colnames(df)[4] <- 'r'
colnames(df)[5] <- 'clidelta'

p <- ggplot(df,aes(x=t))
#
# need to investigate mapping= designator. somehow it's necessary to overlay y-axis coordinated graph.
#

p <- p + geom_point(mapping=aes(y=i,colour=clidelta),stat="identity", position="identity",size=0.8)
p <- p + geom_path(aes(y=i),stat="identity", position="identity",colour="black",linetype="dotted")
p <- p + scale_x_date(date_breaks = "2 year", date_labels = "%Y")
# same as above about mapping=
p <- p + geom_path(aes(y=g),colour='red')
p <- p + geom_path(aes(y=e),colour='blue')
p <- p+annotate("text",label=as.character("10%"),x=as.Date("2020-01-01"), y=550,colour='black')
p <- p + geom_hline(yintercept = 250,size=0.5,linetype=1,colour="white",alpha=1)
p <- p+annotate("text",label=as.character("5%"),x=as.Date("2020-01-01"), y=300,colour='black')
p <- p + geom_hline(yintercept = -250,size=0.5,linetype=1,colour="white",alpha=1)
p <- p+annotate("text",label=as.character("-5%"),x=as.Date("2020-01-01"), y=-200,colour='black')
p <- p + geom_bar(aes(y=r,fill=clidelta),stat = "identity",colour="black") # need identity to draw value itself.

p <- p + theme(axis.title.x=element_blank(),axis.title.y=element_blank())
p <- p + labs(title = "SPX + Theory + Residual% + CLI 1 month Delta",fill="CLI Delta",colour = "CLI Delta")
p <- p +scale_color_brewer(palette="Spectral",na.value = "black",name = "CLI Delta", labels = c("High","mid High","mid Low","Low","NA"))
p <- p +scale_fill_brewer(palette="Spectral",na.value = "black",name = "CLI Delta", labels = c("High","mid High","mid Low","Low","NA"))
p <- p+theme( rect = element_rect(fill = "grey88",
                                  colour = "black",
                                  size = 0,
                                  linetype = 1),
              panel.background = element_rect(fill = "grey88",
                                              colour = "lightblue"),
              axis.title.x=element_blank(),
              axis.title.y=element_blank()
              )
# p <- p + scale_fill_discrete(name="Experimental\nCondition")
# p <- p +scale_colour_discrete(name  ="CLI Delta") + scale_shape_discrete(name  ="CLI Delta")

plot(p)
# remove unnecessary function.
remove(func)
remove(df)



ggplot() histogram by 8 colors.



  • categorize cli delta into 8 groups as to equalize the number of each group as much as possible.
    • count total number.
    • calculate "should be" number of each group.
    • find watermark according to the expected number of iteration.
    • in "func()" and mapply, to put mark "a" to "g" on each month
  • make histogram with  "stack" parameter.
  • remove unnecessary objects.

func <- function(z,a,b,c,d,e,f,g){
  w <- watermark
  x <- z
  if(is.na(x)){return(NA)}
  if(x > a){return("a")}
  if(x > b){return("b")}
  if(x > c){return("c")}
  if(x > d){return("d")}
  if(x > e){return("e")}
  if(x > f){return("f")}
  if(x >= g){return("g")}
  if(x < g){return("h")}
  return(x)
}


# spx_mean <- apply.monthly(SP5[,4],mean)["1970::"]
cov <- apply.monthly(SP5[,4],sd)/apply.monthly(SP5[,4],mean)["1970::"]
lag <- 5
delta <- append(as.vector(diff(cli_xts$oecd,lag)[paste(substr(as.Date(head(index(cov),1)),1,7),"::",sep="")]),rep(NA,length(index(cov)) - length(diff(cli_xts$oecd,lag)[paste(substr(as.Date(head(index(cov),1)),1,7),"::",sep="")])))
# calculate threshold. output 7 thresholds for 8 groups. store into watermark
# use mapply shown in the sample when data.frame() is done.
# when number of colors =8 floor(length(na.omit(delta))/##9##)*seq(2,##8##,1)
watermark <- sort(delta,decreasing = T)[floor(length(na.omit(delta))/9)*seq(2,8,1)]
# mapply(func,delta,watermark[1],watermark[2],watermark[3],watermark[4],watermark[5] )
df <- data.frame(
  i=as.vector(cov),
  # when number of colors =8 watermark[1] till watermark[7] are used.
  c=as.vector(mapply(func,delta,watermark[1],watermark[2],watermark[3],watermark[4],watermark[5],watermark[6],watermark[7])),
  t=as.Date(index(cov)))



p <- ggplot(df, aes(x=i,fill=factor(c)))
# p <- p + geom_histogram(bins=80,position = "stack", alpha = 0.9)
# p <- p + geom_histogram(bins=160,position = "fill", alpha = 0.8)
p <- p + geom_histogram(bins=160,position = "stack", alpha = 0.8)
# p <- p +scale_fill_brewer(na.value = "grey50",name = "CLI Delta", labels = c(as.character(round(watermark,digits=2)),"Less than above","NA"))
legendlable <- c(paste("more than ",as.character(round(watermark,digits=2)),sep=""),"Less than above","NA")
legendlable
p <- p +scale_fill_brewer(palette="Spectral",na.value = "grey50",name = "CLI Delta", labels = legendlable)
p <- p + theme(axis.title.x=element_blank(),axis.title.y=element_blank())
p <- p + theme(panel.background = element_rect(fill = "black",
                                               colour = "lightblue"))
               # panel.grid = element_blank())
# p <- p +scale_fill_brewer(palette="Spectral",na.value = "grey50",name = "CLI Delta", labels = c(as.character(round(watermark,digits=2)),"Less than above","NA"))

plot(p)
# remove(s)
# remove(df)
remove(watermark)
remove(delta)
remove(cov)





2019年9月16日月曜日

ggplot() sample



dataframe を準備

  • 異なる種類のデータを統合するときはこのdata frame作成時にやる。
  • カラム名がggplot()内部での変数名となるので、注意!
df <- data.frame(
  i=as.vector(cov),
  c=as.vector(mapply(func,delta,watermark[1],watermark[2],watermark[3],watermark[4],watermark[5],watermark[6],watermark[7])),
  t=as.Date(index(cov)))

ggplot作成の第一歩

データフレームおよび基礎データの指定

  • データとして使用するデータフレームを指定する。
  • aes()でいろいろな見た目パラメータを指定する。ここではfillで塗りつぶし色を決定する基準を指定している。
  • 離散量で指定したいときはfactor()を使用する。
  • その中で順序を指定したい時はlevels パラメターを使用する。そこでoriginalが取りうる値とその相互の順序を指定する。
    • xlables <- factor(original,levels=c("2-F--64", "2-F-65-", "2-M--64", "2-M-65-", "2-U--64", "2-U-65-", "2-U-UNK","1-F--64", "1-F-65-", "1-M--64", "1-M-65-", "1-M-UNK", "1-U--64", "1-U-65-", "1-U-UNK")) 
  • 連続量を使うときはnumericなどの型のデータをそのまま指定すれば良い。
  • y軸と違いx軸は大概の場合、常に共通なので、aes(x=<column name of data frame>) を使ってここで指定すると良い。
  • 以降はオブジェクト p に対して更新を行う。
    p <- ggplot(df, aes(x=i,fill=factor(c)))

グラフの種類およびyデータの指定

  • ぞれぞれ棒グラフの作成、分散図の作成、線グラフの作成する。
  • aes(y=<column name of data fram>,,,)の形でデータフレームのカラム名のうち適当なものをy軸のデータとして指定する。
  • 色を付けたい場合はaes(fill=<column name of data frame>)などとする。
  • geom_segment()は任意の2つの座標の間に線分を引く
p <- p + geom_bar(aes(y=r),stat = "identity",fill='pink',colour="black") # need identity to draw value itself.
p <- p + geom_point(aes(y=i),stat="identity", position="identity",colour="green",size=0.8)
p <- p + geom_path(aes(y=i),stat="identity", position="identity",colour="black",linetype="dotted")
p <- p + geom_segment(x=as.Date("1985-01-01"),y=log(168),xend=as.Date("2019-09-01"),yend=log(3000),color='white',size=0.02,linetype=2)

棒グラフに色をつけたい場合。

この例ではfillとcolor(colour)双方をつかっている。詳細はscale_color_brewer()scale_fill_brewer()で追加の設定を行う。
p <- p + geom_bar(aes(y=r,fill=clidelta),stat = "identity",colour="black") # need identity to draw value itself.
p <- p + geom_point(mapping=aes(y=i,colour=clidelta),stat="identity", position="identity",size=0.8)
色指定したいデータについて使用するパレットとそれぞれの色値に対する説明を入力する。fillの場合はこのscale_fill_brewer()を使う。 scale_fill_discrete()とかscale_colour_discrete()あるがよくわからないことがまだ多いので、注意して使うこと。
legendlable <- c(paste("more than ",as.character(round(watermark,digits=2)),sep=""),"Less than above","NA")
p <- p +scale_fill_brewer(palette="Spectral",na.value = "grey50",name = "CLI Delta", labels = legendlable)

ヒストグラムの色々

  • 目的に応じて、"identity","stack","fill" それぞれのパラメータを指定する。
  • ヒストグラムのx軸に連続量を指定する場合、離散量を指定する場合で結果が全く異なるので注意すること。
p <- p + geom_histogram(aes(fill=sign),position = "identity", alpha = 0.3,bins=120)
p <- p + geom_histogram(bins=50,position = "fill", alpha = 0.9)
p <- p + geom_histogram(bins=80,position = "stack", alpha = 0.9)

X軸の間隔指定

Date型のときに使用する。
p <- p + scale_x_date(date_breaks = "2 year", date_labels = "%Y")
時間量に応じてグラデーションで色を変えたいときに使用する。
p <- p + scale_fill_date(low = "green3" , high = "darkgreen")


グラフを重ねて書きたい場合

最初にNULLでggplotオブジェクトを作成しそこに順次各種グラフを追加していく。

p <- ggplot(NULL)

df <- data.frame(

  t=last(index(tokyo_death),length_graph),

  value=last(tokyo_death[,1],length_graph)

)

p <- g+geom_bar(data=df, aes(x = t, y = value),stat = "identity",alpha=0.5,colour="red",fill="red")

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)

)

p <- g+geom_line(data=df, aes(x = t, y = value))

plot(g)


タイトルおよび凡例

全体タイトルとレジェンドのタイトルを指定する。共通なタイトルを指定するとレジェンドを一つに統合できる。一方、パレットを使用しないスケールの場合はhueを指定する。

p <- p + labs(title = "SPX + Theory + Residual + CLI Delta",fill="CLI Delta",colour = "CLI Delta")
p <- p + scale_color_brewer(palette="Spectral",na.value = "black",name = "CLI Delta", labels = c("High","mid High","mid Low","Low","NA"))
p <- p + scale_fill_brewer(palette="Spectral",na.value = "black",name = "CLI Delta", labels = c("High","mid High","mid Low","Low","NA"))

p <- p + scale_fill_hue(name='regions')

任意の線を引く

水平線を引く。
p <- p + geom_vline(xintercept=seq(as.Date(paste(substr(index(head(spx_mean,1)),1,7),"-01",sep="")),as.Date("2019-01-01"),by='years'), colour="white",size=0.4,alpha=0.5)
垂直線を引く。
p <- p + geom_vline(xintercept=seq(as.Date(paste(substr(index(head(spx_mean,1)),1,7),"-01",sep="")),as.Date("2019-01-01"),by='years'), colour="white",size=0.4,alpha=0.5)
任意の直線を引く
p <- p + geom_segment(x=as.Date("1985-01-01"),y=log(168),xend=as.Date("2019-09-01"),yend=log(3000),color='white',size=0.02,linetype=2)
回帰線
p <- p + stat_smooth(aes(x=t,y=i),method="loess",color='white',size=0.3)

テーマ

x-y軸のタイトル消去
p <- p + theme(axis.title.x=element_blank(),axis.title.y=element_blank())
台紙の色指定
p <- p + theme(rect = element_rect(fill = "grey88",
                                  colour = "black",
                                  size = 0,
                                  linetype = 1))
パネルの色指定およびグリッドの消去
p <- p + theme(panel.background = element_rect(fill = "grey88",
                                              colour = "lightblue"),
             panel.grid = element_blank())

その他

コメント入力
p <- p + annotate("text",label=as.character(s),x=as.Date("2000-01-01"), y=log(s*1.03),colour='white')
未検証
p <- p + theme(axis.text = element_text(colour = "red", size = rel(1.5)))


パレットの検証
library(RColorBrewer)
display.brewer.all()

過去に使用した関数

  • p <- p + geom_bar
  • p <- p + geom_path
  • p <- p + labs
  • p <- p + scale_fill_date
  • p <- p + scale_x_date
  • p <- p + geom_bar
  • p <- p + geom_histogram
  • p <- p + geom_hline
  • p <- p + geom_path
  • p <- p + geom_point
  • p <- p + geom_vline
  • p <- p + labs
  • p <- p + scale_color_brewer
  • p <- p + scale_fill_brewer
  • p <- p + scale_x_date
  • p <- p + theme
  • p <- p + scale_fill_hue