Draw histgrams to compare high and low

Write function to draw the histgram being optimized for recording blood pressure.

my_bp_hist_x <- function(par_xts,kikan,highlow,len,loc_x,loc_y,br,ymax,xmin,xmax,color){
  # kikan like kikan or "2018::06-21::2018-07-13"
  # highlow 1 or 2. 1 means high, 2 for low.(loc_y+0.04)
  # len is a number of samples.
  # loc_x is x-axis position kikan 0 to 1.
  # loc_y is y-axis.
  # br is breaks for hist()
  # ymax is max value of y-axis.
  # xmin and xmax are for x-axis
  # color is from 1 to 9?
  # function(bp.bangkok,"::2016-06-21,2",length(as.vector(bp.bangkok["2018-06-21::"][,2]),0.3,0.5,20,55,120,2)
  hist(as.vector(last(par_xts[kikan][,highlow],n=len)),breaks=br,xlim=c(xmin,xmax),ylim=c(0,ymax),col=color)
  # axis(side=2, pos=84,labels=F)
  axis(side=2, pos=round(mean(par_xts[kikan][,highlow])),labels=F)
  graph_dim <- par('usr')
  text( graph_dim[1] + graph_dim[2] *  loc_x   ,(graph_dim[4] - graph_dim[3]) * (loc_y+0) + graph_dim[3] ,paste("#",len,sep="="),adj=c(0,0))
  text( graph_dim[1] + graph_dim[2] * loc_x   ,(graph_dim[4] - graph_dim[3]) * (loc_y+0.1) + graph_dim[3] ,paste("mean",round(mean(par_xts[kikan][,highlow]),2),sep="="),adj=c(0,0))
  text( graph_dim[1] + graph_dim[2] * loc_x   ,(graph_dim[4] - graph_dim[3]) * (loc_y+0.2) + graph_dim[3] ,paste("sd",round(sd(par_xts[kikan][,highlow]),2),sep="="),adj=c(0,0))
  text( graph_dim[1] + graph_dim[2] * loc_x   ,(graph_dim[4] - graph_dim[3]) * (loc_y+0.3) + graph_dim[3] ,paste("period ",kikan,sep="="),adj=c(0,0))
}


Use as below which replaces Aug 2 entry.


par(mfrow=c(3,1))
my_bp_hist_x(bp.bangkok,"::2018-06-19",2,70,0.1,0.6,20,15,55,100,4)
my_bp_hist_x(bp.bangkok,"2018-06-21::2018-07-13",2,length(bp.bangkok["2018-06-21::2018-07-13"][,2]),0.1,0.6,20,15,55,100,5)
my_bp_hist_x(bp.bangkok,"2018-07-15::",2,length(bp.bangkok["2018-07-15::"][,2]),0.1,0.6,10,15,55,100,6)
par(mfrow=c(1,1))


par(mfrow=c(3,1))
my_bp_hist_x(bp.bangkok,"::2018-06-19",1,70,0.1,0.6,30,15,90,150,4)
my_bp_hist_x(bp.bangkok,"2018-06-21::2018-07-13",1,length(bp.bangkok["2018-06-21::2018-07-13"][,2]),0.1,0.6,20,15,90,150,5)
my_bp_hist_x(bp.bangkok,"2018-07-15::",1,length(bp.bangkok["2018-07-15::"][,2]),0.1,0.6,20,15,90,150,6)
par(mfrow=c(1,1))

Draw histgram in multiple panes.

1. Draw histgrams in 3 panes. Each pane represents
1. before 2018/06/19
2. between 06/21 and 07/13
3. after 07/15
2. Set the mean value mark by axis command.
3. Output overlay text for Standard Deviation, Mean value and the number of samples.

par(mfrow=c(3,1))
hist(as.vector(last(bp.bangkok["::2018-06-19"][,2],n=length(as.vector(bp.bangkok["2018-06-21::"][,2])))),breaks=20,xlim=c(55,100),ylim=c(0,15),col=3)
# axis(side=2, pos=84,labels=F) 
axis(side=2, pos=round(mean(bp.bangkok["::2018-06-19"][,2])),labels=F) 
graph_dim <- par('usr') 
text( (graph_dim[1] + graph_dim[2]) / 2.2   ,(graph_dim[4] - graph_dim[3]) * 0.64 + graph_dim[3] ,paste("#",round(length(bp.bangkok["2018-06-21::"][,2]),2),sep="="),adj=c(0,0))
text( (graph_dim[1] + graph_dim[2]) / 2.2   ,(graph_dim[4] - graph_dim[3]) * 0.70 + graph_dim[3] ,paste("mean",round(mean(bp.bangkok["::2018-06-19"][,2]),2),sep="="),adj=c(0,0))
text( (graph_dim[1] + graph_dim[2]) / 2.2   ,(graph_dim[4] - graph_dim[3]) * 0.76 + graph_dim[3] ,paste("sd",round(sd(bp.bangkok["::2018-06-19"][,2]),2),sep="="),adj=c(0,0))
hist(as.vector(bp.bangkok["2018-06-21::2018-07-13"][,2]),breaks=20,xlim=c(55,100),ylim=c(0,15),col=2)
axis(side=2, pos=round(mean(bp.bangkok["2018-06-21::2018-07-13"][,2])),labels=F) 
graph_dim <- par('usr') 
text( (graph_dim[1] + graph_dim[2]) / 2.2   ,(graph_dim[4] - graph_dim[3]) * 0.64 + graph_dim[3] ,paste("#",round(length(bp.bangkok["2018-06-21::2018-07-13"][,2]),2),sep="="),adj=c(0,0))
text( (graph_dim[1] + graph_dim[2]) / 2.2   ,(graph_dim[4] - graph_dim[3]) * 0.70 + graph_dim[3] ,paste("mean",round(mean(bp.bangkok["2018-06-21::2018-07-13"][,2]),2),sep="="),adj=c(0,0))
text( (graph_dim[1] + graph_dim[2]) / 2.2   ,(graph_dim[4] - graph_dim[3]) * 0.76 + graph_dim[3] ,paste("sd",round(sd(bp.bangkok["2018-06-21::2018-07-13"][,2]),2),sep="="),adj=c(0,0))
hist(as.vector(bp.bangkok["2018-07-15::"][,2]),breaks=10,xlim=c(55,100),ylim=c(0,15),col=4)
axis(side=2, pos=round(mean(bp.bangkok["2018-07-15::"][,2])),labels=F) 
graph_dim <- par('usr') 
text( (graph_dim[1] + graph_dim[2]) / 1.8   ,(graph_dim[4] - graph_dim[3]) * 0.64 + graph_dim[3] ,paste("#",round(length(bp.bangkok["2018-07-15::"][,2]),2),sep="="),adj=c(0,0))
text( (graph_dim[1] + graph_dim[2]) / 1.8   ,(graph_dim[4] - graph_dim[3]) * 0.70 + graph_dim[3] ,paste("mean",round(mean(bp.bangkok["2018-07-15::"][,2]),2),sep="="),adj=c(0,0))
text( (graph_dim[1] + graph_dim[2]) / 1.8   ,(graph_dim[4] - graph_dim[3]) * 0.76 + graph_dim[3] ,paste("sd",round(sd(bp.bangkok["2018-07-15::"][,2]),2),sep="="),adj=c(0,0))
par(mfrow=c(1,1))

PAYEMS monthly delta and S&P500

plot(to.monthly(SP5)[kikan][,1],axes=T)
par(new=T)
plot(diff(PA[kikan]),axes=F,ylim=c(-900,900),col=3)
par(new=T)
plot(rep(0,length(PA[kikan])),col=2,ylim=c(-900,900),axes=F,type='l')

Calculate weekly standard deviation

Calculate weekly standard deviation from the start of 2017.

> w<-c();for(i in seq(0,trunc(length(day_xts$inc["2017"])/7)-1 ,1)){w<-append(w,sd(day_xts$inc["2017"][i*7 + c(1,2,3,4,5,6,7)]))};print(w);plot(w)

 [1] 65.57148 15.34834 29.44486 12.59251 21.51522 17.22816 10.75042 31.44534 25.34054 27.51017 24.53375 28.72779 31.99702 28.58238 22.88168 28.34818 17.07546
[18] 21.53181 42.79352 32.35444 19.64446 25.16895 24.29580 16.39832

for the case of 2016. The year is already finished, then magic number 52-1 is in place. it is ugly but acceptable for me.

> w<-c();for(i in seq(0,51 ,1)){w<-append(w,sd(day_xts$inc["2016"][i*7 + c(1,2,3,4,5,6,7)]))};print(w)

 [1] 48.606780 25.084619 24.538599 11.443443 20.049938 12.266874 20.221394 14.648257 17.941705 17.201883 19.345234 13.740798 30.581663 15.650955 18.997494 35.491783
[17] 19.474036 24.333986 19.550910 28.180540 18.927179 16.340135 14.404034 21.351926 23.831752 20.868066 22.401105 23.999008 20.497387 27.022918 28.678846 22.003247
[33] 29.384965 25.602269 17.305380 25.746937 20.846377 27.993196 22.112268 26.525819 24.013885 29.664794 22.166040 33.642377 21.883893 36.787032 24.006943 26.567076
[49] 20.317715  7.653197 14.893271 21.084863

2 graphs S&P500 theoretical vs. actual and PAYEMS moving average

when kg is as below
> kg
[1] "1992-01-01::2017-03-31"

kikan <- kg
plot(predict(lm(to.quarterly(SP5[kikan])[,1] ~ apply.quarterly(PA[kikan],mean) * apply.quarterly(UC[kikan],mean) * GDP[kikan] - apply.quarterly(UC[kikan],mean) )),type='l',col=2,axes=F,ylim=c(0,2500),ylab="")
par(new=T)
plot(to.quarterly(SP5[kikan])[,1],type='l',ylim=c(0,2500))

l <- 300  # length of sample period
a <- 6     # base for moving average
plot(last(seq(as.Date("1992-01-01"),as.Date("2017-05-01"),by="months"),n=(l-a)),na.trim(filter(diff(last(PA,n=l)),rep(1,a))/a),type='h',ylim=c(-800,300))

Below will overlay moving average upon actual. Don't forget to set 'l' and 'a'.

 plot(last(seq(as.Date("1992-01-01"),as.Date("2017-05-01"),by="months"),n=(l-a)),na.trim(filter(diff(last(PA,n=l)),rep(1,a))/a),type='l',ylim=c(-800,500))
 par(new=T)
 plot(last(diff(PA),n=(l-a)),ylim=c(-800,500),col=2,type='h')

check normal distribution (shapilo.test)

Outlook often gives viewers a false impression. Shapiro.test function will test whether data are according to a normal distribution. When p-Value is more than 0.05, the distribution can be called "normal".

Samples below shows when data are extracted just from Monday to Friday without considering national holiday, it's not a normal distribution. However, other groups 1) workday (business day only), 2) all-day and 3) weekend (Saturday and Sunday) distribution are normal. The degree of each distribution is different, though.


> shapiro.test(i[weekday])

Shapiro-Wilk normality test
data:  i[weekday]
W = 0.97156, p-value = 0.02788

> shapiro.test(i[workday])

Shapiro-Wilk normality test
data:  i[workday]
W = 0.98811, p-value = 0.5692

> shapiro.test(i)

Shapiro-Wilk normality test

data:  i

W = 0.98489, p-value = 0.1182

> shapiro.test(i[weekend])

Shapiro-Wilk normality test

data:  i[weekend]

W = 0.95236, p-value = 0.07851