S&P500 and GDP and other macro-economical data - Case-Shiller 10-City Composite Home Price Index

getSymbols("GDP",src='FRED')
getSymbols("PAYEMS",src="FRED")
getSymbols("UNDCONTSA",src="FRED")

getSymbols('SPCS10RSA',src='FRED')
CS <- SPCS10RSA

k2k <- "2000-01-01::2017-12-01"

PAq <- apply.quarterly(PA[k2k],mean)
UCq <- apply.quarterly(UC[k2k],mean)
CSq <- apply.quarterly(CS[k2k],mean)
G <- GDP

summary(lm(apply.quarterly(SP5[k2k],mean)[,1] ~ PAq[k2k] * UCq[k2k] * G[k2k]*CSq[k2k] - UCq[k2k] -G[k2k] - PAq[k2k]*G[k2k] - UCq[k2k]*G[k2k]*CSq[k2k]))

> summary(lm(apply.quarterly(SP5[k2k],mean)[,1] ~ PAq[k2k] * UCq[k2k] * G[k2k]*CSq[k2k] - UCq[k2k] -G[k2k] - PAq[k2k]*G[k2k] - UCq[k2k]*G[k2k]*CSq[k2k]))

Call:

lm(formula = apply.quarterly(SP5[k2k], mean)[, 1] ~ PAq[k2k] * 

    UCq[k2k] * G[k2k] * CSq[k2k] - UCq[k2k] - G[k2k] - PAq[k2k] * 

    G[k2k] - UCq[k2k] * G[k2k] * CSq[k2k])

Residuals:

     Min       1Q   Median       3Q      Max 

-193.182  -33.564   -0.483   37.984  186.693 

Coefficients:

                                    Estimate Std. Error t value Pr(>|t|)    

(Intercept)                       -1.080e+03  3.794e+02  -2.846  0.00591 ** 

PAq[k2k]:UCq[k2k]                  7.692e-05  6.148e-06  12.512  < 2e-16 ***

PAq[k2k]:CSq[k2k]                 -2.019e-04  3.278e-05  -6.160 5.13e-08 ***

PAq[k2k]:UCq[k2k]:G[k2k]          -5.157e-09  4.380e-10 -11.776  < 2e-16 ***

PAq[k2k]:UCq[k2k]:CSq[k2k]        -1.639e-07  3.181e-08  -5.153 2.59e-06 ***

PAq[k2k]:G[k2k]:CSq[k2k]           2.161e-08  1.608e-09  13.442  < 2e-16 ***

PAq[k2k]:UCq[k2k]:G[k2k]:CSq[k2k]  9.053e-12  1.822e-12   4.968 5.19e-06 ***

---

Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 75.55 on 65 degrees of freedom

Multiple R-squared:  0.9713, Adjusted R-squared:  0.9686 

F-statistic: 366.6 on 6 and 65 DF,  p-value: < 2.2e-16

> my_sp5cs(k2k,19981.13, 147810,1120,221)

[1] "m_m params! apply.quarter - UC w/ nominal GDP"

[1] 2725.425

k = [1] "2000-01-01::2017-12-01"

plot(xts(as.numeric(residuals(lm(apply.quarterly(SP5[k2k],mean)[,1] ~ PAq[k2k] * UCq[k2k] * G[k2k]*CSq[k2k] - UCq[k2k] -G[k2k] - PAq[k2k]*G[k2k] - UCq[k2k]*G[k2k]*CSq[k2k]))),seq(as.Date(strsplit(k2k,"::")[[1]][1]),as.Date(strsplit(k2k,"::")[[1]][2]),by="quarters")),type='p')

Plot actual and theoretical - S&P500

plot(merge(apply.quarterly(SP5[k2k],mean)[,1],predict(lm(apply.quarterly(SP5[k2k],mean)[,1] ~ PAq[k2k] * UCq[k2k] * G[k2k]*CSq[k2k] - UCq[k2k] -G[k2k] - PAq[k2k]*G[k2k] - UCq[k2k]*G[k2k]*CSq[k2k]))))

plot(xts(residuals(lm(apply.quarterly(SP5[k2k],mean)[,1] ~ PAq[k2k] * UCq[k2k] * G[k2k]*CSq[k2k] - UCq[k2k] -G[k2k] - PAq[k2k]*G[k2k] - UCq[k2k]*G[k2k]*CSq[k2k]))[,1],seq(as.Date(strsplit(k2k,"::")[[1]][1]),as.Date(strsplit(k2k,"::")[[1]][2]),by='quarters')),type='p')

New S&P500 = GDP and other macro-economical data

getSymbols("GDP",src='FRED')
getSymbols("PAYEMS",src="FRED")
getSymbols("UNDCONTSA",src="FRED")
getSymbols('SPCS20RSA',src='FRED')

CS <- SPCS20RSA
k2k <- "2000-01-01::2017-12-01"
PAq <- apply.quarterly(PA[kikan],mean)
UCq <- apply.quarterly(UC[kikan],mean)
CSq <- apply.quarterly(CS[kikan],mean)

summary(lm(apply.quarterly(SP5[k2k],mean)[,1] ~ PAq[k2k] * UCq[k2k] * G[k2k]*CSq[k2k] - UCq[k2k] -G[k2k] - PAq[k2k]*G[k2k] - UCq[k2k]*G[k2k]*CSq[k2k]))


Residuals:
     Min       1Q   Median       3Q      Max
-208.815  -29.865    1.545   39.597  203.971

Coefficients:
                                    Estimate Std. Error t value Pr(>|t|) 
(Intercept)                       -1.250e+03  4.070e+02  -3.071  0.00311 **
PAq[k2k]:UCq[k2k]                  8.127e-05  6.593e-06  12.327  < 2e-16 ***
PAq[k2k]:CSq[k2k]                 -2.300e-04  3.765e-05  -6.109 6.27e-08 ***
PAq[k2k]:UCq[k2k]:G[k2k]          -5.465e-09  4.813e-10 -11.355  < 2e-16 ***
PAq[k2k]:UCq[k2k]:CSq[k2k]        -1.731e-07  3.413e-08  -5.070 3.54e-06 ***
PAq[k2k]:G[k2k]:CSq[k2k]           2.503e-08  1.927e-09  12.988  < 2e-16 ***
PAq[k2k]:UCq[k2k]:G[k2k]:CSq[k2k]  9.177e-12  1.954e-12   4.697 1.41e-05 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 80.57 on 65 degrees of freedom
Multiple R-squared:  0.9674, Adjusted R-squared:  0.9643
F-statistic:   321 on 6 and 65 DF,  p-value: < 2.2e-16



my_sp5cs(k2k,last(G[k2k]),as.vector(last(PAq[k2k])),as.vector(last(UCq[k2k])),as.vector(last(CSq[k2k])))

S&P500,GDP,PAYEMS,UNDCONTSA 2018JAN version

# set parameters.
kg <- "1992-01-01::2017-12-31"
kikan <- "1992-01-01::2017-12-31"
l <- 48  # # of months to predict
r <- 1.05 # pesumed GDP growth rate
# download data from FRE.
getSymbols("GDP",src="FRED")
getSymbols("PAYEMS",src="FRED")
getSymbols("UNDCONTSA",src="FRED")
PA <- PAYEMS
UC <- UNDCONTSA
SP5 <- as.xts(read.zoo(read.csv("~/SP5.csv")))
G <- GDP
i <- seq(2,l/3,1)  # seq of quarters to predict
d <- as.Date(as.yearqtr(seq(Sys.Date(),as.Date("2100-12-31"),by="quarters")[i])) # pick up the first day of each quarters.
m_GDP <- as.xts(as.vector(last(GDP)) * r**(i/4),d)
m_PA <- (as.xts(forecast(auto.arima(PA),h=l)$mean[1:l],as.Date(as.yearmon(seq(mondate(index(last(PA)))+1,by=1,length.out=l))))[(3-month(index(last(PA))) %% 3) + seq(1,l-3,3)])[d]
m_UC <- (as.xts(forecast(auto.arima(UC),h=l)$mean[1:l],as.Date(as.yearmon(seq(mondate(index(last(UC)))+1,by=1,length.out=l))))[(3-month(index(last(UC))) %% 3) + seq(1,l-3,3)])[d]
my_sp5(kg,m_GDP,m_PA,m_UC)

> kikan
[1] "1992-01-01::2017-12-31"
> summary(lm(apply.quarterly(SP5[kikan],mean)[,1] ~ apply.quarterly(PA[kikan],mean) * apply.quarterly(UC[kikan],mean) * G[kikan] - apply.quarterly(UC[kikan],mean) -G[kikan] ))

Call:
lm(formula = apply.quarterly(SP5[kikan], mean)[, 1] ~ apply.quarterly(PA[kikan],
    mean) * apply.quarterly(UC[kikan], mean) * G[kikan] - apply.quarterly(UC[kikan],
    mean) - G[kikan])

Residuals:
    Min      1Q  Median      3Q     Max
-255.30  -74.08  -11.65   71.56  300.72


                                                                           Pr(>|t|) 
(Intercept)                                                                2.53e-05 ***
apply.quarterly(PA[kikan], mean)                                           1.13e-05 ***
apply.quarterly(PA[kikan], mean):apply.quarterly(UC[kikan], mean)          4.89e-10 ***
apply.quarterly(PA[kikan], mean):G[kikan]                                  3.12e-09 ***
apply.quarterly(UC[kikan], mean):G[kikan]                                   < 2e-16 ***
apply.quarterly(PA[kikan], mean):apply.quarterly(UC[kikan], mean):G[kikan]  < 2e-16 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 117.4 on 98 degrees of freedom
Multiple R-squared:  0.9505, Adjusted R-squared:  0.9479
F-statistic: 376.1 on 5 and 98 DF,  p-value: < 2.2e-16

draw the graph to record blood pressure

# draw the graph data are picked up when above 101 to avoid anormaly.
#
plot(bp.xts[,-3][bp.xts$High > 101],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 > 101]))),index(bp.xts[bp.xts$High > 101])),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 > 101]))),index(bp.xts[bp.xts$High > 101])),ylim=c(60,160),on=1,col=5,lwd=1)
# add aux axis
axis(2,at=c(135,130,85))


#
# prepare xts data
#
# read data as csv format and convert to xts to plot
# bp <- read.csv("~/bp.csv")
bp <- read.csv("~/Downloads/bp - シート1.csv")
#
# remove the CSV file after the update
# don't use "~"  $HOME is okay R's "system" command is very awkward
# on this front.
#
system("rm \"$HOME/Downloads/bp - シート1.csv\"")
#
# below is also possible of course.
# system("rm \"/Users/honomoto/Downloads/bp - シート1.csv\"")
# while below is okay
# system("rm ~/Downloads/bp*.csv")
# but ~ is not usable in "" as comma suppress meta character.
# system("rm \"~/Downloads/bp - シート1.csv\"")

# create xts, remove data and time from the original and
# append POSIX time stamp.
bp.xts <- xts(bp[,c(-1,-2)],as.POSIXct(paste(bp$Date,bp$Time,sep=" ")))
# calculate weekly mean except the data of the sample whose $High is <=100
#
apply.weekly(bp.xts[bp.xts$High > 101],mean)

manage blood pressure part.2

# read data as csv format and convert to xts to plot
# bp <- read.csv("~/bp.csv")
bp <- read.csv("~/Downloads/bp - シート1.csv")
#

# create xts, remove data and time from the original and
# append POSIX time stamp.
bp.xts <- xts(bp[,c(-1,-2)],as.POSIXct(paste(bp$Date,bp$Time,sep=" ")))
# calculate weekly mean except the data of the sample whose $High is <=100
#
apply.weekly(bp.xts[bp.xts$High > 101],mean)
# drow plot with standard line at 85 and 130.

plot(merge(bp.xts,cbind(rep(85,length(bp.xts[,1])),rep(130,length(bp.xts[,1])))),col = c("red", "blue","green","green"),lwd=c(3,3,2,2),major.ticks='days',grid.ticks.on='days')

# or exclude the data $High is <= 100

plot(merge(bp.xts[,-3][bp.xts$High > 101],cbind(rep(85,length(bp.xts[,1][bp.xts$High > 101])),rep(130,length(bp.xts[,1][bp.xts$High > 101])))),col = c("red", "blue","green","green"),lwd=c(3,3,2,2),major.ticks='days',grid.ticks.on='days')


plot(xts(bp[,c(-1,-2,-5)],as.POSIXct(paste(bp$Date,bp$Time,sep=" "))))



# calculate 4 days moving average

 plot(last(index(xts(bp[,c(-1,-2,-5)],as.POSIXct(paste(bp$Date,bp$Time,sep=" ")))),n=length(na.omit(filter(bp$High,rep(1,4))/4))),na.omit(filter(bp$High,rep(1,4))/4),type='l')

# same but for 7 days.

 plot(last(index(xts(bp[,c(-1,-2,-5)],as.POSIXct(paste(bp$Date,bp$Time,sep=" ")))),n=length(na.omit(filter(bp$High,rep(1,7))/7))),na.omit(filter(bp$High,rep(1,7))/7),type='l')


# exclude the data whose High is more than 150 and less than 101 and calculate
# mean value.
mean(na.omit(bp.xts[bp.xts$High < 150 & bp.xts$High > 101][,1]))

# 1) draw High and Low by red and blue each
# 2) define linewidth for High, Low and rep(130,74)
# 3) draw an additional line at 130
# 4) set x lim tick frequency to "days"

plot(merge(bp.xts,rep(130,74)),col = c("red", "blue","black"),lwd=c(3,3,0.5),major.ticks='days')

# 6) and add grid frequency to "days"

plot(merge(bp.xts,rep(130,74)),col = c("red", "blue","green"),lwd=c(3,3,0.5),major.ticks='days',grid.ticks.on='days')

plot(merge(bp.xts,cbind(rep(85,length(bp.xts[,1])),rep(130,length(bp.xts[,1])))),col = c("red", "blue","green","green"),lwd=c(3,3,2,2),major.ticks='days',grid.ticks.on='days')