2021年5月10日月曜日

横並び棒グラフ dodge tidyr gather 三指数

 
> last(dmdf,3)
    01Hokkaido 02Aomori 03Iwate 04Miyagi 05Akita 06Yamagata 07Fukushima 08Ibaraki 09Tochigi 10Gunma 11Saitama 12Chiba 13Tokyo 14Kanagawa
477          4        0       0        1       0          0           1         0         1       2         1       7       6          6
478          3        2       1        0       0          0           0         0         0       0         1       5       6          2
479          7        0       1        2       0          1           1         0         0       0         2       1       3          4
    15Niigata 16Toyama 17Ishikawa 18Fukui 19Yamanashi 20Nagano 21Gifu 22Shizuoka 23Aichi 24Mie 25Shiga 26Kyoto 27Osaka 28Hyogo 29Nara
477         0        1          1       0           0        2      1          1       3     1       3       0      50      39      4
478         0        0          4       0           0        0      0          0       2     0       0       2      41       6      0
479         0        0          1       0           0        0      0          0       5     0       0       0      19       8      0
    30Wakayama 31Tottori 32Shimane 33Okayama 34Hiroshima 35Yamaguchi 36Tokushima 37Kagawa 38Ehime 39Kochi 40Fukuoka 41Saga 42Nagasaki
477          0         0         0         3           0           0           1        0       3       0         4      0          0
478          1         0         0         4           1           0           2        0       1       0         1      0          1
479          0         0         0         2           0           2           1        0       1       0         3      0          0
    43Kumamoto 44Oita 45Miyazaki 46Kagoshima 47Okinawa          t
477          1      0          0           0         1 2021-05-07
478          0      0          0           0         0 2021-05-08
479          0      0          0           0         0 2021-05-09

> head(w)
    tokyo osaka hyogo          t
400    11     2     1 2021-02-19
401    27     4     2 2021-02-20
402    17     1     0 2021-02-21
403     9     3     1 2021-02-22
404    11     5     3 2021-02-23
405    17     4     7 2021-02-24

> dplyr::filter(df,t > as.Date("2021-05-07"))
           t variable value
1 2021-05-08    tokyo     6
2 2021-05-09    tokyo     3
3 2021-05-08    osaka    41
4 2021-05-09    osaka    19
5 2021-05-08    hyogo     6
6 2021-05-09    hyogo     8



w <- data.frame(tokyo=dmdf[,13],osaka=dmdf[,27],hyogo=dmdf[,28],t=dmdf$t)
w <- last(w,80)
df <-w  %>% tidyr::gather(variable,value,-t)  #exclude t to preserve time stamp or use "df <-w  %>% tidyr::gather(variable,value,1:3)"
# df
p <- ggplot(df,aes(x=t, y=value, color=variable,fill=variable)) 
# p <- p+geom_histogram(stat = "identity",position = "stack")
p <- p+geom_bar(stat="identity",width=1,position = "dodge")
plot(p)
# end here

# if "df <-w  %>% tidyr::gather(variable,value)  #exclude t to preserve time stamp"
# 警告メッセージ: 
# attributes are not identical across measure variables;
# they will be dropped 
# see w[1:3,] contents are numeric  and w$t is Date.



w <- data.frame(tokyo=dmdf[,13],osaka=dmdf[,27],hyogo=dmdf[,28],t=dmdf$t)
w <- last(w,80)
df <-w  %>% tidyr::gather(variable,value,-t) #exclude t to preserve time stamp or use "df <-w  %>% tidyr::gather(variable,value,1:3)"
# df
p <- ggplot(df,aes(x=t, y=value,fill=variable)) 
p <- p +scale_fill_brewer(palette="Accent")
p <- p+geom_bar(stat="identity",width=1,position = "dodge")
p <- p + theme(panel.background = element_rect(fill = "grey60",
                                               colour = "lightblue"))
plot(p)









len <- 20
w <- data.frame(spx=as.vector(last(weeklyReturn(GSPC),len)),ndx=as.vector(last(weeklyReturn(NDX),len)),dji= as.vector(last(weeklyReturn(DJI),len))    ,t=index(last(weeklyReturn(GSPC),len)))
df <-w  %>% tidyr::gather(index,weeklyreturn,-t)
df$t[len*seq(1,3,1)] <- df$t[dim(df)[1]-1]+7  # adjust each index's last entry. 3 came from gspc, ndx and di.
for(i in seq(1,len*3,1)){if((length(seq(as.Date("1970-01-05"),as.Date(df$t[i]),by='days'))%%7) != 5){df$t[i] <- df$t[i]+5-length(seq(as.Date("1970-01-05"),as.Date(df$t[i]),by='days'))%%7 }} # for the case national holiday distort date.
p <- ggplot(df,aes(x=t, y=weeklyreturn, color=index,fill=index))
p <- p + xlab("") + ylab("週間収益率")
# p <- p + geom_point(aes(y=case_per_capita,size=sign,color=r),alpha=1)
# p <- p+annotate("text",label=pref_db[,1],x=df[,2], y=df[,1]+0.1,colour='black',family = "HiraKakuProN-W3",size=3)
p <- p + theme_gray (base_family = "HiraKakuPro-W3")
p <- p + scale_color_brewer(name="指数",labels=c("ダウ30種","ナスダック","S&P500"),palette="Accent")
p <- p + scale_fill_brewer(name="指数",labels=c("ダウ30種","ナスダック","S&P500"),palette="Accent")
p <- p+geom_bar(stat="identity",width=7,position = "dodge")
plot(p)





2021年4月29日木曜日

名目GDP四半期伸び率 vs. CLI Delta Nominal GDP quater basis legend title #spline #gdp

 

名目GDP四半期伸び率 vs. CLI Delta Nominal GDP quater basis を対照プロットする。


period <- "2001::2021-06"  # end period month should be march,june,september or december.
w <- (GDPC1/lag(GDPC1))[period] - 1
result <- lm(as.vector(diff(apply.quarterly(cli_xts$oecd,mean))[period]) ~ w)
df <- data.frame(delta=as.vector(diff(apply.quarterly(cli_xts$oecd,mean))[period]) ,gdp=as.vector(w)
                 ,spx= as.vector(quarterlyReturn(GSPC)[period])  )

p <- ggplot(df, aes(y=delta,x=gdp))
p <- p + theme_gray (base_family = "HiraKakuPro-W3")
p <- p + xlab("GDP四半期伸び率") + ylab("CLI delta")
p <- p + geom_point(alpha=1,aes(color=spx))
p <- p + scale_color_gradient(low = "red", high = "green")
p <- p + geom_smooth(method = "lm",se=F,size=0.5)
p <- p + geom_vline(xintercept =(-1*result$coefficients[1] / result$coefficients[2]),
                    size=0.5,linetype=2,colour="red",alpha=0.5)
p <- p + geom_hline(yintercept = last(as.vector(diff(apply.quarterly(cli_xts$oecd,mean))[period])),
                    size=0.5,linetype=2,colour="green",alpha=0.5)
p <- p + geom_vline(xintercept = last(w),
                    size=0.5,linetype=2,colour="green",alpha=0.5)
p <- p+annotate("text",label=as.character(round(-1*result$coefficients[1] / result$coefficients[2],4)),x=(-1*result$coefficients[1] / result$coefficients[2]), y=-5.1,colour='red')
p <- p+annotate("text",label=as.character(round(last(w),4)),x=last(w), y=-5.1,colour='green')

p <- p +  labs(color="S&P500\n四半期\n収益率")
plot(p)
summary(result)

spline を使用して補完の上月次GDPを算出その上でグラフを作成する。

period <- "2015::2021-09"  # end period month should be march,june,september or december.

w <- spline(GDPC1[period],xmin=1/3,xmax=length(GDPC1[period])+1/3,n=3*length(GDPC1[period])+1,method='natural')$y
# w <- w/lag(w) -1
w <- w[2:length(w)]/w[1:(length(w)-1)] - 1
# w <- (GDPC1/lag(GDPC1))[period] - 1

# result <- lm(as.vector(diff(cli_xts$oecd)[period]) ~ w)
result <- lm(as.vector(diff(cli_xts$oecd)[period]) ~ w)
df <- data.frame(delta=as.vector(diff((cli_xts$oecd))[period]) ,gdp=as.vector(w)
                 ,spx= as.vector(monthlyReturn(GSPC)[period])  )

p <- ggplot(df, aes(y=delta,x=gdp))
p <- p + theme_gray (base_family = "HiraKakuPro-W3")
p <- p + xlab("GDP補正月間伸び率") + ylab("CLI delta")
p <- p + geom_point(alpha=1,aes(color=spx))
p <- p + scale_color_gradient(low = "red", high = "green")
p <- p + geom_smooth(method = "lm",se=F,size=0.5)
p <- p + geom_vline(xintercept =(-1*result$coefficients[1] / result$coefficients[2]),
                    size=0.5,linetype=2,colour="red",alpha=0.5)
p <- p + geom_hline(yintercept = last(as.vector(diff((cli_xts$oecd))[period])),
                    size=0.5,linetype=2,colour="green",alpha=0.5)
p <- p + geom_vline(xintercept = last(w),
                    size=0.5,linetype=2,colour="green",alpha=0.5)
p <- p+annotate("text",label=as.character(round(-1*result$coefficients[1] / result$coefficients[2],4)),x=(-1*result$coefficients[1] / result$coefficients[2]), y=-5.1,colour='red')
p <- p+annotate("text",label=as.character(round(last(w),4)),x=last(w), y=-5.1,colour='green')

p <- p +  labs(color="S&P500\n月間\n収益率")
plot(p)
summary(result)



> summary(result)

Call:
lm(formula = as.vector(diff(apply.quarterly(cli_xts$oecd, mean))[period]) ~ w)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.44728 -0.33122  0.01372  0.19903  2.09477 

Coefficients:
            Estimate Std. Error t value Pr(>|t|)    
(Intercept) -0.45260    0.06811  -6.645 3.53e-09 ***
w           47.61859    3.71852  12.806  < 2e-16 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 0.5276 on 79 degrees of freedom
Multiple R-squared:  0.6749, Adjusted R-squared:  0.6708 
F-statistic:   164 on 1 and 79 DF,  p-value: < 2.2e-16

2021年4月27日火曜日

EPS 2021APR27

 


MacBook-Pro-16:~/Downloads$ cat eps.txt 
12/31/2022 $54.64 $50.47 20.39 21.88 $202.75 $188.97
9/30/2022 $52.43 $49.07 21.08 22.63 $196.14 $182.75
6/30/2022 $48.86 $45.76 21.84 23.43 $189.30 $176.48
3/31/2022 $46.82 $43.67 22.64 24.32 $182.63 $170.06
12/31/2021 $48.03 $44.25 23.22 24.79 $178.09 $166.77
9/30/2021 $45.59 $42.80 24.58 26.86   $168.24 $153.96
6/30/2021 $42.19 $39.33 25.76 28.69   $160.55 $144.14
3/31/2021 (25.4%) 3972.89 $42.28 $40.39 28.49 33.72 $145.15 $122.64
12/31/2020 3756.07 $38.18 $31.44 30.69 39.90   $122.37 $94.13
9/30/2020 3363.00 $37.90 $32.98 27.26 34.24   $123.37 $98.22
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


MacBook-Pro-16:~/Downloads$  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] <- 98.22
eps_year_xts["2019::"][8] <- 94.13
eps_year_xts["2019::"][9] <- 122.64
eps_year_xts["2019::"][10] <- 144.14
eps_year_xts["2019::"][11] <- 153.96
eps_year_xts["2019::"][12] <- 166.77
eps_year_xts["2019::"][13] <- 170.06
eps_year_xts["2019::"][14] <- 176.48
eps_year_xts["2019::"][15] <- 182.75
eps_year_xts["2019::"][16] <- 188.97

> eps_year_xts["2020::"]
             [,1]
2020-01-01 116.33
2020-04-01  99.23
2020-07-01  98.22
2020-10-01  91.15
2021-01-01 114.36
2021-04-01 133.35
2021-07-01 141.08
2021-10-01 155.56

> eps_year_xts <- append(eps_year_xts,as.xts(rep(0,4),seq(as.Date("2022-01-01"),as.Date("2022-10-01"),by='quarters')))


2021年4月24日土曜日

Create your own discrete scale Source: R/scale-manual.r, R/zxx.r

original is here.


Create your own discrete scaleSource: R/scale-manual.r, R/zxx.r

These functions allow you to specify your own set of mappings from levels in the data to aesthetic values.

  • scale_colour_manual(..., values, aesthetics = "colour", breaks = waiver())
  • scale_fill_manual(..., values, aesthetics = "fill", breaks = waiver())
  • scale_size_manual(..., values, breaks = waiver())
  • scale_shape_manual(..., values, breaks = waiver())
  • scale_linetype_manual(..., values, breaks = waiver())
  • scale_alpha_manual(..., values, breaks = waiver())
  • scale_discrete_manual(aesthetics, ..., values, breaks = waiver())


Arguments passed on to discrete_scale

  • palette A palette function that when called with a single integer argument (the number of levels in the scale) returns the values that they should take (e.g., scales::hue_pal()).
    • limits One of:
    • NULL to use the default scale values
    • A character vector that defines possible values of the scale and their order
    • A function that accepts the existing (automatic) values and returns new ones
  • drop Should unused factor levels be omitted from the scale? The default, TRUE, uses the levels that appear in the data; FALSE uses all the levels in the factor.
  • na.translate  Unlike continuous scales, discrete scales can easily show missing values, and do so by default. If you want to remove missing values from a discrete scale, specify na.translate = FALSE.
  • na.value  If na.translate = TRUE, what aesthetic value should the missing values be displayed as? Does not apply to position scales where NA is always placed at the far right.
  • scale_name  The name of the scale that should be used for error messages associated with this scale.
  • name The name of the scale. Used as the axis or legend title. If waiver(), the default, the name of the scale is taken from the first mapping used for that aesthetic. If NULL, the legend title will be omitted.
  • labels 
    • One of:
    • NULL for no labels
    • waiver() for the default labels computed by the transformation object
    • A character vector giving labels (must be same length as breaks)
    • A function that takes the breaks as input and returns labels as output
  • guide A function used to create a guide or its name. See guides() for more information.
  • super The super class to use for the constructed scale
  • values a set of aesthetic values to map data values to. The values will be matched in order (usually alphabetical) with the limits of the scale, or with breaks if provided. If this is a named vector, then the values will be matched based on the names instead. Data values that don't match will be given na.value.
  • aesthetics Character string or vector of character strings listing the name(s) of the aesthetic(s) that this scale works with. This can be useful, for example, to apply colour settings to the colour and fill aesthetics at the same time, via aesthetics = c("colour", "fill").
  • breaks
    • One of:
    • NULL for no breaks
    • waiver() for the default breaks (the scale limits)
    • A character vector of breaks
    • A function that takes the limits as input and returns breaks as output


Details


The functions scale_colour_manual(), scale_fill_manual(), scale_size_manual(), etc. work on the aesthetics specified in the scale name: colour, fill, size, etc. However, the functions scale_colour_manual() and scale_fill_manual() also have an optional aesthetics argument that can be used to define both colour and fill aesthetic mappings via a single function call (see examples). The function scale_discrete_manual() is a generic scale that can work with any aesthetic or set of aesthetics provided via the aesthetics argument.


Color Blindness

Many color palettes derived from RGB combinations (like the "rainbow" color palette) are not suitable to support all viewers, especially those with color vision deficiencies. Using viridis type, which is perceptually uniform in both colour and black-and-white display is an easy option to ensure good perceptive properties of your visulizations. 

The colorspace package offers functionalitiesto generate color palettes with good perceptive properties,to analyse a given color palette, like emulating color blindness,and to modify a given color palette for better perceptivity.

Examples

p <- ggplot(mtcars, aes(mpg, wt)) +
  geom_point(aes(colour = factor(cyl)))
p + scale_colour_manual(values = c("red", "blue", "green"))

# It's recommended to use a named vector
cols <- c("8" = "red", "4" = "blue", "6" = "darkgreen", "10" = "orange")
p + scale_colour_manual(values = cols)

# You can set color and fill aesthetics at the same time
ggplot(
  mtcars,
  aes(mpg, wt, colour = factor(cyl), fill = factor(cyl))
) +
  geom_point(shape = 21, alpha = 0.5, size = 2) +
  scale_colour_manual(
    values = cols,
    aesthetics = c("colour", "fill")
  )

# As with other scales you can use breaks to control the appearance
# of the legend.
p + scale_colour_manual(values = cols)
p + scale_colour_manual(
  values = cols,
  breaks = c("4", "6", "8"),
  labels = c("four", "six", "eight")
)

# And limits to control the possible values of the scale
p + scale_colour_manual(values = cols, limits = c("4", "8"))
#> Warning: Removed 7 rows containing missing values (geom_point).
p + scale_colour_manual(values = cols, limits = c("4", "6", "8", "10"))


2021年4月23日金曜日

palette color 色の定義 scale_fill_manual manal scale マニュアルスケール 色見本


1. define gradient private palette and pick up 5 colors. 

my_palette <- colorRampPalette(c("#FF0000","#FFFF00","#00FF00","#00FFFF","#0000FF"))
plot_col <- my_palette(5)[5:1]


2.pick up 5 colors from the standard palette "Spectral".

RColorBrewer::brewer.pal(5,"Spectral")[5:1]

3.display sample

barplot(rep(1,15), col=my_palette(15), axes=FALSE)

barplot(rep(1,11), col=RColorBrewer::brewer.pal(11,"Spectral")[11:1], axes=FALSE)

*note that spectral max color num is 11.

3. define private palette

my_palette <- colorRampPalette(c("#FF0000","#FFFF00","#00FF00","#00FFFF","#0000FF"))
# create 47 colors from my_palette.
g <- g + scale_fill_manual(name='regions',values=my_palette(47),labels=pref_jp)
# can't use scale_fill_hue with 'values=' parameter.



2021年4月11日日曜日

atan plot.default pch

 
> atan2(diff(cli_xts$oecd),diff(diff(cli_xts$oecd))) %>% last(.,)
               oecd
2021-03-01 1.495797

adjust to the date of the last available data.

v <- atan2(diff(cli_xts$oecd),diff(diff(cli_xts$oecd))) %>% last(.,240)
w <- monthlyReturn(GSPC)[
"::2021-03"] %>% last(.,240) # update if necessary.
par(bg="grey",fg="white")
plot.default(v,w,type='p')
tmp <- par('usr')
plot.new()
plot.default(v,w,xlim=c( tmp[1],tmp[2]), ylim=c(tmp[3], tmp[4]),type='p')
abline(v=pi/2)
abline(v=0)
par(new=T)
i <- 12
plot.default(last(v,i),last(w,i),xlim=c( tmp[1],tmp[2]), ylim=c(tmp[3], tmp[4]),type='b',pch='+',col='brown')
for(k in seq(1,9,1)){
  j <- k %% 8 + 1
  par(new=T)
  print(i)
  print(j)
  plot.default(last(v,i)[k],last(w,i)[k],xlim=c( tmp[1],tmp[2]), ylim=c(tmp[3], tmp[4]),type='p',pch=as.character(k),col=j)
}



2021年4月3日土曜日

tidyverse tidyr dplyr detach search searchpath library library(help= )

 


> search()   
> searchpaths()    # インストール場所の絶対パスを得る

# あるライブラリ中のオブジェクトの一覧を見る library(help=パッケージ名) 
# もちろんそのライブラリが既にインストールされている必要があります。

> library(gstat)       # ライブラリ gstat をロード
> library(help=gstat) # gstat 中のオブジェクト(関数、データセット)の一覧を得る


> library(tidyverse)
─ Attaching packages ─────────────────────────────────────────────────── tidyverse 1.3.0 ─
✓ tibble  3.1.0     ✓ dplyr   1.0.5
✓ tidyr   1.1.3     ✓ stringr 1.4.0
✓ readr   1.4.0     ✓ forcats 0.5.1
✓ purrr   0.3.4     
─ Conflicts ───────────────────────────────────────────────────── tidyverse_conflicts() ─
x stringr::boundary() masks strucchange::boundary()
x dplyr::filter()     masks stats::filter()
x dplyr::first()      masks xts::first()

x dplyr::lag()        masks stats::lag()
x dplyr::last()       masks xts::last()
x dplyr::select()     masks MASS::select()

> data(mtcars)
> mtcars <- as_tibble(mtcars, rownames = "model") %>% mutate(cyl = as.character(cyl))
> g <- ggplot(mtcars, aes(x = mpg, y = wt, color = cyl)) +
+     geom_text(aes(label = model), family = family_sans, fontface = "plain") +
+     labs(x =paste(family_serif, "ボールドを使用"), y = paste(family_serif, "イタリックを使用"),
+          title = paste(family_serif, "ボールドイタリックを使用")) +
+     annotate("text", x = 10, y = 2, label = paste(family_sans, "標準書体を使用"), hjust = 0) +
+     theme(
+         text = element_text(family = family_serif, face = "plain"),
+         title = element_text(face = "bold.italic"),
+         axis.title = element_text(face = "italic"),
+         axis.title.x = element_text(face = "bold")
+     )
> g
Error in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y,  : 
  polygon edge not found
In addition: There were 50 or more warnings (use warnings() to see the first 50)
> g <- ggplot(mtcars, aes(x = mpg, y = wt, color = cyl)) +
+     geom_text(aes(label = model), family = family_sans, fontface = "plain") +
+     labs(x =paste(family_serif, "ボールドを使用"), y = paste(family_serif, "イタリックを使用"),
+          title = paste(family_serif, "ボールドイタリックを使用")) +
+     annotate("text", x = 10, y = 2, label = paste(family_sans, "標準書体を使用"), hjust = 0) +
+     theme(
+         text = element_text(family = family_serif, face = "plain"),
+         title = element_text(face = "bold.italic"),
+         axis.title = element_text(face = "italic"),
+         axis.title.x = element_text(face = "bold")
+     )
> g
Error in grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y,  : 
  polygon edge not found
In addition: Warning message:
In grid.Call(C_textBounds, as.graphicsAnnot(x$label), x$x, x$y,  :
  no font could be found for family "Hiragino Mincho ProN"
>
last(dmdf[,13],10)   # blocked by dplyr
Error in order(order_by)[[n]] : subscript out of bounds
> detach("package:tidyverse", unload=TRUE)
> last(dmdf[,13],10)
Error in order(order_by)[[n]] : subscript out of bounds
> search()
 [1] ".GlobalEnv"           "package:forcats"      "package:stringr"      "package:dplyr"        "package:purrr"        "package:readr"       
 [7] "package:tidyr"        "package:tibble"       "package:RColorBrewer" "package:ggplot2"      "tools:rstudio"        "package:datasets"    
[13] "package:beepr"        "package:forecast"     "package:mondate"      "package:vars"         "package:lmtest"       "package:urca"        
[19] "package:strucchange"  "package:sandwich"     "package:MASS"         "package:utils"        "package:graphics"     "package:grDevices"   
[25] "package:quantmod"     "package:TTR"          "package:xts"          "package:zoo"          "package:stats"        "package:methods"     
[31] "Autoloads"            "package:base"        
> detach("package:dplyr", unload=TRUE)
Warning message:
‘dplyr’ namespace cannot be unloaded:
  namespace ‘dplyr’ is imported by ‘dbplyr’, ‘broom’, ‘tidyr’ so cannot be unloaded 
> detach("package:tidyr", unload=TRUE)
Warning message:
‘tidyr’ namespace cannot be unloaded:
  namespace ‘tidyr’ is imported by ‘broom’ so cannot be unloaded 
> detach("package:broom", unload=TRUE)
Error in detach("package:broom", unload = TRUE) : invalid 'name' argument
> detach("package:tidyr", unload=TRUE)
Error in detach("package:tidyr", unload = TRUE) : invalid 'name' argument
> detach("package:dplyr", unload=TRUE)
Error in detach("package:dplyr", unload = TRUE) : invalid 'name' argument
> last(dmdf[,13],10)
 [1]  7 18  6  7 15 16 20 12 10 23

2021年4月1日木曜日

annotate , xlab , ylab , scale_color_hue, size, legend


df <- data.frame(case_per_capita=as.vector(apply(mdf[,-48],2,sum) / pref_db$x2017),pop_density=pref_db$x2017/pref_db$size,sign=pref_db$x2017,r=pref_db[,2])

p <- ggplot(df, aes(x=pop_density,y=case_per_capita,size=sign,color=r))
p <- p +
xlab("人口密度") + ylab("人口あたり件数")
p <- p + geom_point(alpha=1)
p <- p+annotate("text",label=pref_db[1,3],x=df[1,2], y=df[1,1]+0.1,colour='red',family = "HiraKakuProN-W3")
p <- p+annotate("text",label=pref_db[11,3],x=df[11,2], y=df[11,1]+0.1,colour='red',family = "HiraKakuProN-W3")
p <- p+annotate("text",label=pref_db[12,3],x=df[12,2], y=df[12,1]+0.1,colour='red',family = "HiraKakuProN-W3")
p <- p+annotate("text",label=pref_db[13,3],x=df[13,2], y=df[13,1]+0.1,colour='red',family = "HiraKakuProN-W3")
p <- p+annotate("text",label=pref_db[14,3],x=df[14,2], y=df[14,1]+0.1,colour='red',family = "HiraKakuProN-W3")
p <- p+annotate("text",label=pref_db[26,3],x=df[26,2], y=df[26,1]+0.1,colour='red',family = "HiraKakuProN-W3")
p <- p+annotate("text",label=pref_db[23,3],x=df[23,2], y=df[23,1]+0.1,colour='red',family = "HiraKakuProN-W3")
p <- p+annotate("text",label=pref_db[27,3],x=df[27,2], y=df[27,1]+0.1,colour='red',family = "HiraKakuProN-W3")
p <- p+annotate("text",label=pref_db[28,3],x=df[28,2], y=df[28,1]+0.1,colour='red',family = "HiraKakuProN-W3")
p <- p+annotate("text",label=pref_db[40,3],x=df[40,2], y=df[40,1]+0.1,colour='red',family = "HiraKakuProN-W3")
p <- p+annotate("text",label=pref_db[47,3],x=df[47,2], y=df[47,1]+0.1,colour='red',family = "HiraKakuProN-W3")
# p <- p +scale_color_brewer(palette="Spectral")
# p <- p +scale_color_brewer(palette=rainbow(47))
# p <- p +scale_color_brewer()
p <- p + theme_gray (base_family = "HiraKakuPro-W3")
p <- p + scale_color_hue(name="都道府県",labels=pref_db[,3])
# p <- p + guides(fill = guide_legend(reverse = F,order = 2),label = TRUE)
p <- p + guides(size = guide_legend(title="人口"))
# don't forget to set "color=". otherwise fails to show up.
p <- p + geom_smooth(method = "lm",se=F,color="red",size=1)
# p + scale_colour_manual(values = pref_db[,3])
plot(p)



for the final version, see below.


palette("Alpabet")
df <- data.frame(case_per_capita=as.vector(apply(mdf[,-48],2,sum) / pref_db$x2017),pop_density=pref_db$x2017/pref_db$size,sign=pref_db$x2017,r=pref_db[,2])
# p <- ggplot(df, aes(x=pop_density,y=case_per_capita,size=sign,color=r))
p <- ggplot(df, aes(x=pop_density,y=case_per_capita,size=sign,color=r))
p <- p + xlab("人口密度") + ylab("人口あたり件数")
p <- p + geom_point(alpha=1)
p <- p+annotate("text",label=pref_db[attributes(df[df$case_per_capita > 3 ,])$row.names,3],x=df[attributes(df[df$case_per_capita > 3 ,])$row.names,2], y=df[attributes(df[df$case_per_capita > 3 ,])$row.names,1]+0.1,colour='red',family = "HiraKakuProN-W3")
# p <- p +scale_color_brewer(palette="Spectral")
# p <- p +scale_color_brewer(palette=rainbow(47))
# p <- p +scale_color_brewer()
p <- p + theme_gray (base_family = "HiraKakuPro-W3")
p <- p + scale_color_hue(name="都道府県",labels=pref_db[,3])
# p <- p + guides(fill = guide_legend(reverse = F,order = 2),label = TRUE)
p <- p + guides(size = guide_legend(title="人口"))
# don't forget to set "color=". otherwise fails to show up.
p <- p + geom_smooth(method = "lm",se=F,color="red",size=1)  
plot(p)

df <- data.frame(death_per_capita=as.vector(apply(dmdf[,-48],2,sum) / pref_db$x2017),pop_density=pref_db$x2017/pref_db$size,sign=pref_db$x2017,r=pref_db[,2])
# p <- ggplot(df, aes(x=pop_density,y=case_per_capita,size=sign,color=r))
p <- ggplot(df, aes(x=pop_density,y=death_per_capita,size=sign,color=r))
p <- p + xlab("人口密度") + ylab("人口あたり死者数")
p <- p + geom_point(alpha=1)
p <- p+annotate("text",label=pref_db[attributes(df[df$death_per_capita > 0.07 ,])$row.names,3],x=df[attributes(df[df$death_per_capita > 0.07 ,])$row.names,2], y=df[attributes(df[df$death_per_capita > 0.07 ,])$row.names,1]+0.002,colour='red',family = "HiraKakuProN-W3")

# p <- p +scale_color_brewer(palette="Spectral")
# p <- p +scale_color_brewer(palette=rainbow(47))
# p <- p +scale_color_brewer()
p <- p + theme_gray (base_family = "HiraKakuPro-W3")
p <- p + scale_color_hue(name="都道府県",labels=pref_db[,3])
# p <- p + guides(fill = guide_legend(reverse = F,order = 2),label = TRUE)
p <- p + guides(size = guide_legend(title="人口"))
# don't forget to set "color=". otherwise fails to show up.
p <- p + geom_smooth(method = "lm",se=F,color="red",size=1)  
plot(p)

2021年3月22日月曜日

upgrading to R 4

 # library(magrittr) ??
install.packages("tidyverse")
install.packages("tidyr")
install.packages("beepr")
install.packages("forecast")
install.packages("mondate")
install.packages("vars")
install.packages("quantmod")
install.packages("xts")

upgrade version up

2021年3月14日日曜日

金利 NASDAQ 相関 cor.test() intersect()


choose date both TNX and NDX delta is available.


 d <- intersect(index(diff(TNX[,4])[!is.na(diff(TNX[,4]))]), index(diff(NDX[,4])[!is.na(diff(NDX[,4]))]))

check outputs.


head(as.Date(d))

[1] "2007-01-04" "2007-01-05" "2007-01-09" "2007-01-10" "2007-01-11" "2007-01-12"

Check correlation in 2021.


cor.test(diff(TNX[,4])["2021"][as.Date(d)],diff(NDX[,4])["2021"][as.Date(d)])

Pearson's product-moment correlation
data:  diff(TNX[, 4])["2021"][as.Date(d)] and diff(NDX[, 4])["2021"][as.Date(d)]
t = -2.3358, df = 36, p-value = 0.02519
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 -0.61153962 -0.04874962
sample estimates:
       cor 
-0.3627793 
 警告メッセージ: 
 z + c(-1, 1) * sigma * qnorm((1 + conf.level)/2) で: 
  Recycling array of length 1 in array-vector arithmetic is deprecated.
  Use c() or as.vector() instead.

Check correlation all time.


cor.test(diff(TNX[,4])[as.Date(d)],diff(NDX[,4])[as.Date(d)])

Pearson's product-moment correlation
data:  diff(TNX[, 4])[as.Date(d)] and diff(NDX[, 4])[as.Date(d)]
t = 9.819, df = 2787, p-value < 2.2e-16
alternative hypothesis: true correlation is not equal to 0
95 percent confidence interval:
 0.1467388 0.2184915
sample estimates:
      cor 
0.1828587 
 警告メッセージ: 
 z + c(-1, 1) * sigma * qnorm((1 + conf.level)/2) で: 
  Recycling array of length 1 in array-vector arithmetic is deprecated.
  Use c() or as.vector() instead.


> plot.default(diff(TNX[,4])[as.Date(d)],diff(NDX[,4])[as.Date(d)])
> plot.default(diff(TNX[,4])["2021"][as.Date(d)],diff(NDX[,4])["2021"][as.Date(d)])

2021年3月13日土曜日

NASDAQ EPS

download data from here.

modify data like below. use   %s/\s/,/g   in vi.

2020-12-31,1.34
2020-09-30,1.58
2020-06-30,1.45
2020-03-31,1.22
2019-12-31,1.21
       <skip>
2005-06-30,0.13
2005-03-31,0.13

store data into csv file. " ~/nasdaqeps.txt" is used in the sample. change end date accordingly!.

w <- as.xts(read.zoo(read.csv("~/nasdaqeps.txt",header=FALSE)))
ndx_eps <- as.xts(as.vector(na.omit(filter(as.vector(w),rep(1,4)))),seq(as.Date("2005-10-01"),as.Date("2020-10-01"),by='quarters'))

when k2007 is as below.

>  k2007
[1] "2007-01-01::2020-12-31"

summary(lm(apply.quarterly(NDX[,4][k2007],mean) ~ ndx_eps[k2007]+apply.quarterly(PA[k2007],mean)+apply.quarterly(UC[k2007],mean)+apply.quarterly(CS[k2007],mean)))

Call:
lm(formula = apply.quarterly(NDX[, 4][k2007], mean) ~ ndx_eps[k2007] + 
    apply.quarterly(PA[k2007], mean) + apply.quarterly(UC[k2007], 
    mean) + apply.quarterly(CS[k2007], mean))

Residuals:
   Min     1Q Median     3Q    Max 
-787.4 -293.1 -105.6  106.1 2281.4 

Coefficients:
                                   Estimate Std. Error t value Pr(>|t|)    
(Intercept)                      -1.335e+04  8.966e+02 -14.893  < 2e-16 ***
ndx_eps[k2007]                    2.612e+02  7.864e+01   3.321 0.001664 ** 
apply.quarterly(PA[k2007], mean)  7.181e-01  5.575e-02  12.882  < 2e-16 ***
apply.quarterly(UC[k2007], mean) -2.917e+00  9.409e-01  -3.100 0.003144 ** 
apply.quarterly(CS[k2007], mean)  4.479e+01  1.096e+01   4.088 0.000154 ***
---
Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1

Residual standard error: 532.2 on 51 degrees of freedom
Multiple R-squared:  0.9611, Adjusted R-squared:  0.958 
F-statistic: 314.8 on 4 and 51 DF,  p-value: < 2.2e-16


2021年3月10日水曜日

gawk awk 正規表現 エスケープ sed gsed

 正規表現メモ

gawk 正規表現


\w
これは単語を構成する任意のキャラクタ、つまり 文字、数字、それとアンダースコアにマッチする演算子である。 これは [[:alnum:]_] の簡潔な表現とみなして良い。
\W
これは単語を構成する要素にならない任意のキャラクタにマッチする 演算子である。これは [^[:alnum:]_] の簡潔な表現とみなして良い。
\<
これは単語の先頭にある空文字列にマッチする演算子である。 例えば、/\<away/`away'にマッチするが、 `stowaway'にはマッチしない。
\>
これは単語の末尾にある空文字列にマッチする演算子である。 例えば、/stow\>/`stow'にマッチするが、 `stowaway'にはマッチしない。
\y
これは単語の先頭、あるいは末尾の空文字列とマッチする演算子である (つまり語の区切りとマッチするということである)。 例えば、`\yballs?\y'は独立した単語として `ball' にも `balls'にもマッチする。注意! 他のGNU softwareは単語の区切りに\bを使うらしい。
\B
この演算子は単語中の空文字列にマッチする。言い換えると、`\B'は二つ の単語の構成要素文字の間にある空文字列にマッチするということである。例え ば、/\Brat\B/ は`crate'にマッチする。しかし、`dirty rat' にはマッチしない。`B'は簡単にいうと`\y'の反対語である。
\( .... \) () でマッチの塊を扱う
\n  \n は n 番目の () に対応 ([a-z]*) \1*

gsed 正規表現


\< 
        “単語”の先頭。
\> 
        “単語”の末尾。
\b 
        “単語”の先頭、あるいは末尾。言い換えれば、 “語の区切り”とマッチする。
\B 
        “単語”の中にある空文字列にマッチします
\`  
        バッファの先頭。
\'  
        バッファの末尾。
\w 
        “単語”を構成するキャラクタ、つまり文字と数字それにアンダースコア のいずれか。
\W 
        “単語”を構成するキャラクタ以外のキャラクタの いずれか。
\s 
        「空白」一文字とマッチする。localeの影響に注意。
\S
        「空白以外」の一文字とマッチ。

2021年3月9日火曜日

data.frame locale xts format

 
>  last(PI,8)

                PI
2020-06-01 20032.7
2020-07-01 20173.9
2020-08-01 19624.3
2020-09-01 19762.2
2020-10-01 19627.8
2020-11-01 19386.2
2020-12-01 19499.2
2021-01-01 21453.9

上記がある時日付のフォーマットだけを変更する。

Sys.setlocale("LC_TIME","C") 
data.frame(t=(PI  %>% last(.,8) %>%  index() %>% format(., "%Y %b")), d=(PI  %>% last(.,8) %>%  as.vector()))  
Sys.setlocale("LC_TIME", "ja_JP.UTF-8")


YYYY MMM にするのはこちら

Sys.setlocale("LC_TIME","C")

data.frame(t=(PI  %>% last(.,8) %>%  index() %>% format(., "%b %Y")), d=(PI  %>% last(.,8) %>%  as.vector()))

Sys.setlocale("LC_TIME", "ja_JP.UTF-8") 


         t       d
1 2020 Jun 20032.7
2 2020 Jul 20173.9
3 2020 Aug 19624.3
4 2020 Sep 19762.2
5 2020 Oct 19627.8
6 2020 Nov 19386.2
7 2020 Dec 19499.2
8 2021 Jan 21453.9

>   system("locale")

LANG="ja_JP.UTF-8"
LC_COLLATE="ja_JP.UTF-8"
LC_CTYPE="ja_JP.UTF-8"
LC_MESSAGES="ja_JP.UTF-8"
LC_MONETARY="ja_JP.UTF-8"
LC_NUMERIC="ja_JP.UTF-8"
LC_TIME="ja_JP.UTF-8"
LC_ALL=

日付 フォーマット date time format POSIX

 

1.時刻フィールド


%H     時 (00..23)
%I     時 (01..12)
%k     時 (0..23)
%l     時 (1..12)
%M     分 (00..59)
%p     AM あるいは PM のロケール
%r     時刻、12 時間 (hh:mm:ss [AP]M)
%s     1970-01-01 00:00:00 UTC からの秒数 (標準外の拡張)
%S     秒 (00..60)
%T     時刻、24 時間 (hh:mm:ss)
%X     ロケールによる時刻の表現 (%H:%M:%S)
%Z     タイムゾーン (例 EDT)、あるいはタイムゾーンが決定できないならば無し

2.日付フィールド


%a     ロケールの省略形の曜日名 (Sun..Sat)
%A     ロケールの完全表記の曜日名、可変長 (Sunday..Saturday)
%b     ロケールの省略形の月名 (Jan..Dec)
%B     ロケールの完全表記の月名、可変長 (January..December)
%c     ロケールの日付と時刻 (Sat Nov 04 12:02:33 EST 1989)
%d     月内通算日数 (01..31)
%D     日付 (mm/dd/yy)
%h     %b と同じ
%j     年内通算日数 (001..366)
%m     月 (01..12)
%U     日曜日を週の最初の日とした年内通算週 (00..53)
%w     週のうちの曜日 (0..6) (0 が日曜日)
%W     月曜日を週の最初の日とした年内通算週 (00..53)
%x     ロケールの日付表現 (mm/dd/yy)
%y     年の最後の 2 つの数字 (00..99)
%Y     年 (1970...)