compare PAYEMS and S&P500

Use classic plot function.

• Don't forget to use as.Date to convert index to x-axis, otherwise axis will loose away.
• Locate "axis" sentence in the script with a full caution. It may work at other places but may not either.

plot.default(as.Date(index(to.monthly(SP5)[kikan][,1])),to.monthly(SP5)[kikan][,1],axes=T,type='l',ylim=c(400,3000))
axis(side=2, pos=as.Date("2017-01-01"))
par(new=T)
plot.default(as.Date(index(to.monthly(SP5)[kikan][,1])),rep(2700,length(index(to.monthly(SP5)[kikan][,1]))),type='l',col=3,ylim=c(400,3000))
par(new=T)
plot.default(diff(PA[kikan]),axes=F,ylim=c(-900,900),col=2,type='h')

> to.monthly(SP5)[kikan][200]
        SP5.Open SP5.High SP5.Low SP5.Close  SP5.Volume
 8 2008  1269.42  1313.15 1247.45   1282.83 86266010000

plot.default(as.Date(index(to.monthly(SP5)[kikan][,1])),to.monthly(SP5)[kikan][,1],axes=T,type='l',ylim=c(400,3000),first=grid(28,10),tck=0.01)

par(new=T)

plot.default(diff(PA[kikan]),axes=F,ylim=c(-900,900),col=2,type='h')

abline(v=200)   # draw v-line at 2008 August

par(new=T)

plot.default(as.Date(index(to.monthly(SP5)[kikan][,1])),rep(1282,length(index(to.monthly(SP5)[kikan][,1]))),type='l',col=3,ylim=c(400,3000))

Manage Blood Pressure data

# read data as csv format and convert to xts to plot
#
Sys.setenv(TZ=Sys.timezone())
bp <- read.csv("~/Downloads/bp - シート1.csv")
system("rm \"$HOME/Downloads/bp - シート1.csv\"")
bp.xts <- xts(bp[,c(-1,-2,-6)],as.POSIXct(paste(bp$Date,bp$Time,sep=" "),tz=Sys.timezone()),tz=Sys.timezone())
# weekly average
apply.weekly(bp.xts[bp.xts$High > 95],mean)
#
#  draw the graph
#
plot(bp.xts[,c(1,2)][bp.xts$High > 95],col = c("red", "blue"),lwd=c(3,3,2,2),major.ticks='days',grid.ticks.on='days',type='p',ylim=c(60,160))
# draw a horizontal line at 130 as the benchmark
addSeries(xts(rep(130,length(index(bp.xts[bp.xts$High > 95]))),index(bp.xts[bp.xts$High > 95])),ylim=c(60,160),on=1,col=5,lwd=1)
# draw another line at 85.
addSeries(xts(rep(85,length(index(bp.xts[bp.xts$High > 95]))),index(bp.xts[bp.xts$High > 95])),ylim=c(60,160),on=1,col=5,lwd=1)
# add aux axis
axis(2,at=c(135,130,125,120,115,85))
axis(4,at=c(135,130,125,120,115,85))
# draw the horizontal line at the average
addSeries(xts(rep(mean(bp.xts[,1][bp.xts$High > 95]),length(index(bp.xts[bp.xts$High > 95]))),index(bp.xts[bp.xts$High > 95])),ylim=c(60,160),on=1,col=6,lwd=1)
addSeries(xts(rep(mean(bp.xts[,2][bp.xts$High > 95]),length(index(bp.xts[bp.xts$High > 95]))),index(bp.xts[bp.xts$High > 95])),ylim=c(60,160),on=1,col=6,lwd=1)
events <- xts(c("natrix","weight","abort natrix"),as.Date(c("2018-06-20", "2018-07-14","2018-08-09")))

addEventLines(events, srt=90, pos=2,col=10)

# prepare data according to system timezone. "Asia/Tokyo" in most cases.
#
bp.day <- merge(as.xts(as.vector(bp.xts[,1]),as.Date(index(bp.xts),tz=tzone(bp.xts))),as.vector(bp.xts[,2]))
colnames(bp.day)[1] <- "high"
colnames(bp.day)[2] <- "low"
#
# prepare timezone 2 hours behind "Asia/Tokyo".
#
bp.bangkok <- merge(as.xts(as.vector(bp.xts[,1]),as.Date(index(bp.xts),tz="Asia/Bangkok")),as.vector(bp.xts[,2]))
colnames(bp.bangkok)[1] <- "high"
colnames(bp.bangkok)[2] <- "low"
apply.weekly(bp.bangkok,mean)


#
# t test between before 2018-06-20 and after 2018-06-21
#  2018-06-20 is to add natrix into prescription.
#  2018-07-14 is to start weight training.
#
t.test(as.vector(last(bp.bangkok["::2018-06-20"][,2],n=length(as.vector(bp.bangkok["2018-06-21::"][,2])))),as.vector(bp.bangkok["2018-06-21::"][,2]))
#
# t test no.2
#
t.test(as.vector(bp.bangkok["2018-06-21::2018-07-13"][,2]),as.vector(bp.bangkok["2018-07-15::"][,2]))
#
# t test no.3
#
t.test(as.vector(last(bp.bangkok["::2018-07-13"][,2],n=length(as.vector(bp.bangkok["2018-07-15::"][,2])))),as.vector(bp.bangkok["2018-07-15::"][,2]))
#
# t test no.4
#
t.test(as.vector(last(bp.bangkok["::2018-06-19"][,2],n=length(as.vector(bp.bangkok["2018-08-16::"][,2])))),as.vector(bp.bangkok["2018-08-16::"][,2]))

#
# use anonymous function in line to calculate the average of the data sorted by the time to take samples of higher pressure.
#
mapply(function(x,y){return(mean(na.omit(bp.xts[strptime(format(index(bp.xts),"%H:%M:%S"),"%H:%M:%S") > strptime(x,"%H:%M:%S") & strptime(format(index(bp.xts),"%H:%M:%S"),"%H:%M:%S") < strptime(y,"%H:%M:%S")])
[,1]))},c("05:00:00","06:00:00","07:00:00","08:00:00","09:00:00","10:00:00","11:00:00","23:00:00","00:00:00"),c("05:00:00","06:59:00","07:59:00","08:59:00","09:59:00","10:59:00","11:59:00","23:59:00","00:59:00"))
#
# mean value of the lower prssure.
#
mapply(function(x,y){return(mean(na.omit(bp.xts[strptime(format(index(bp.xts),"%H:%M:%S"),"%H:%M:%S") > strptime(x,"%H:%M:%S") & strptime(format(index(bp.xts),"%H:%M:%S"),"%H:%M:%S") < strptime(y,"%H:%M:%S")])
[,2]))},c("05:00:00","06:00:00","07:00:00","08:00:00","09:00:00","10:00:00","11:00:00","23:00:00","00:00:00"),c("05:00:00","06:59:00","07:59:00","08:59:00","09:59:00","10:59:00","11:59:00","23:59:00","00:59:00"))
#
#  count the number of samples
#
mapply(function(x,y){return(length(na.omit(bp.xts[strptime(format(index(bp.xts),"%H:%M:%S"),"%H:%M:%S") > strptime(x,"%H:%M:%S") & strptime(format(index(bp.xts),"%H:%M:%S"),"%H:%M:%S") < strptime(y,"%H:%M:%S")])
[,1]))},c("05:00:00","06:00:00","07:00:00","08:00:00","09:00:00","10:00:00","11:00:00","23:00:00","00:00:00"),c("05:00:00","06:59:00","07:59:00","08:59:00","09:59:00","10:59:00","11:59:00","23:59:00","00:59:00"))


par(mfrow=c(3,1))
hist(as.vector(last(bp.bangkok["::2018-06-20"][,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-20"][,2])),labels=F) 
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) 
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]) / 2   ,(graph_dim[4] - graph_dim[3]) * 0.70 + graph_dim[3] ,paste("mean",mean(bp.bangkok["2018-07-15::"][,2]),sep="="),adj=c(0,0))
text( (graph_dim[1] + graph_dim[2]) / 2   ,(graph_dim[4] - graph_dim[3]) * 0.75 + graph_dim[3] ,paste("sd",sd(bp.bangkok["2018-07-15::"][,2]),sep="="),adj=c(0,0))
par(mfrow=c(1,1))

IRON ORE VIX and USDJPY

Download monthly Iron Ore price from this URL.

Dont' forget cut&paste table data into the terminal and format them by the editor

cat <INPUT CSV FILE NAME> |sed 's/,//' |awk '{print $1"-"$2","$3","$4}' >  <OUTPUT CSV FILE NAME>

IRON_ORE <- as.xts(as.numeric(read.csv("~/icloud/R/iron_ore_price.csv")[2:240,2]),seq(as.Date("1998-07-01"),as.Date("2018-05-01"),by="months"))

summary(lm(to.monthly(YJUSDJPY)["2008::2017"][,4] ~ to.monthly(VIX)["2008::2017"][,4] + to.monthly(IRON_ORE)["2008::2017"][,4]))

> summary(lm(to.monthly(YJUSDJPY)["2008::2017"][,4] ~ to.monthly(VIX)["2008::2017"][,4] + to.monthly(IRON_ORE)["2008::2017"][,4]))



Call:

lm(formula = to.monthly(YJUSDJPY)["2008::2017"][, 4] ~ to.monthly(VIX)["2008::2017"][, 4] + to.monthly(IRON_ORE)["2008::2017"][, 4])

Residuals:
    Min      1Q  Median      3Q     Max
-32.025  -9.121   0.750   7.764  25.368

Coefficients:
                                                     Estimate Std. Error t value Pr(>|t|) 
(Intercept)                                    98.2969     3.0042  32.720  < 2e-16 ***
to.monthly(VIX)["2008::2017"][, 4]       -0.5002     0.1147  -4.361 2.80e-05 ***
to.monthly(IRON_ORE)["2008::2017"][, 4]   0.1470     0.0238   6.179 9.69e-09 ***

---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1
Residual standard error: 11.57 on 117 degrees of freedom
Multiple R-squared:  0.316, Adjusted R-squared:  0.3043
F-statistic: 27.02 on 2 and 117 DF,  p-value: 2.249e-10



plot(merge(to.monthly(YJUSDJPY["2008::2017"])[,4],predict(lm(to.monthly(YJUSDJPY)["2008::2017"][,4] ~ to.monthly(VIX)["2008::2017"][,4] + to.monthly(IRON_ORE)["2008::2017"][,4]))))

IRON ORE,COPPER,VIX and USDJPY

> summary(lm(to.monthly(YJUSDJPY)["2008::2017"][,4] ~ to.monthly(VIX)["2008::2017"][,4] + to.monthly(COPPER)["2008::2017"][,4] +to.monthly(IRON_ORE)["2008::2017"][,4]))

Call:
lm(formula = to.monthly(YJUSDJPY)["2008::2017"][, 4] ~ to.monthly(VIX)["2008::2017"][,
    4] + to.monthly(COPPER)["2008::2017"][, 4] + to.monthly(IRON_ORE)["2008::2017"][,
    4])

Residuals:
    Min      1Q  Median      3Q     Max
-18.569  -7.629  -0.686   4.995  22.341

Coefficients:
                                          Estimate Std. Error t value Pr(>|t|) 
(Intercept)                              1.426e+02  6.451e+00  22.103  < 2e-16 ***
to.monthly(VIX)["2008::2017"][, 4]      -6.714e-01  9.754e-02  -6.883 3.18e-10 ***
to.monthly(COPPER)["2008::2017"][, 4]   -5.174e-03  6.956e-04  -7.438 1.92e-11 ***
to.monthly(IRON_ORE)["2008::2017"][, 4]  6.416e-02  2.261e-02   2.838  0.00536 **
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 9.565 on 116 degrees of freedom
Multiple R-squared:  0.5369, Adjusted R-squared:  0.5249
F-statistic: 44.82 on 3 and 116 DF,  p-value: < 2.2e-16

plot(merge(to.monthly(YJUSDJPY["2008::2017"])[,4],predict(lm(to.monthly(YJUSDJPY)["2008::2017"][,4] ~ to.monthly(VIX)["2008::2017"][,4] + to.monthly(COPPER)["2008::2017"][,4] +to.monthly(IRON_ORE)["2008::2017"][,4]))))

Calculate your mortgage.

CHECK THIS AS WELL

For the case interest rate is 0.875%, months to pay is 26 full years, and the remaining principle is 50,814,556JPY.

See aDFyear for each year and aDFmonth for each month.

> mortgage(P=50814556, I=0.875, L=(26*12), amort=T, plotData=T)

The payments for this loan are:


Monthly payment: ￥182,154.3 (stored in monthPay)
Total cost: ￥56,832,131

The amortization data for each of the 312 months are stored in "aDFmonth".

The amortization data for each of the 312 years are stored in "aDFyear".

> aDFyear
   Amortization Annual_Payment Annual_Principal Annual_Interest Year
1      50814556        2185851          1748224       437627.33    1
2      49066332        2185851          1763582       422268.88    2
3      47302750        2185851          1779076       406775.49    3
4      45523674        2185851          1794705       391146.00    4

                                     <skip>

25      4332107        2185851          2156580        29270.81   25
26      2175526        2185851          2175526        10324.87   26

Upgrade packages in R

When unable to upgrade packages as they are already loaded, check "package" panel and uncheck loaded ones.

Unchecking the mark is equal to the command "detach("package:RMySQL", unload=TRUE)".
Do same for "DBI".
*Note: Skip "restart R session before upgrade" as it will bring back unloaded packages.

R_proj$cat ~/.Rprofile 

source( "~/R_proj/startup.R" )

R_proj$cat ~/R_proj/startup.R

library(xts)

library(quantmod)

library(vars)

library(mondate)

# library(RMySQL) #junbi

library(forecast)

library(beepr)

Don't forget to comment and uncomment packages in initial start files as required.

Manipulate timezone

#
# prepare data according to system timezone. "Asia/Tokyo" in most cases.
#
bp.day <- merge(as.xts(as.vector(bp.xts[,1]),as.Date(index(bp.xts),tz=tzone(bp.xts))),as.vector(bp.xts[,2]))
colnames(bp.day)[1] <- "high"
colnames(bp.day)[2] <- "low"
#
# prepare timezone 2 hours behind "Asia/Tokyo".
#
bp.bangkok <- merge(as.xts(as.vector(bp.xts[,1]),as.Date(index(bp.xts),tz="Asia/Bangkok")),as.vector(bp.xts[,2]))
colnames(bp.bangkok)[1] <- "high"
colnames(bp.bangkok)[2] <- "low"
apply.weekly(bp.bangkok,mean)

• Set timezone to other than Tokyo, as blood pressure may be measured after the mid night.
• Setting time zone as "Bangkok" allows samples taken from 00:00 to 02:00 to be taken into the day before.

T test to compare before and after the new drug prescribed.

t.test(as.vector(last(bp.bangkok["::2018-06-20"][,2],n=length(as.vector(bp.bangkok["2018-06-21::"][,2])))),as.vector(bp.bangkok["2018-06-21::"][,2]))
par(mfrow=c(2,1))
hist(as.vector(last(bp.bangkok["::2018-06-20"][,2],n=length(as.vector(bp.bangkok["2018-06-21::"][,2])))),breaks=20,xlim=c(55,100),ylim=c(0,10),col=3)
hist(as.vector(bp.bangkok["2018-06-21::"][,2]),breaks=20,xlim=c(55,100),ylim=c(0,10),col=2)
par(mfrow=c(1,1))
#
# mean value before
#
mean(last(bp.bangkok["::2018-06-20"][,2],n=length(as.vector(bp.bangkok["2018-06-21::"][,2]))))
#
# mean value after
#
mean(as.vector(bp.bangkok["2018-06-21::"][,2]))

• 2018-06-20 is the day to start to add the new drug for me.
• Comparing low blood pressure before and after the date.
• The number of sample is standardized.
• Drawing histograms to compare graphically.