merge all data eps spline gpuc

merge all data

result.gpuc <- 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])
result.eps <- lm(apply.quarterly(SP5[,4][k2k],mean) ~ eps_year_xts[k2k]+apply.quarterly(PA[k2k],mean)+apply.quarterly(CS[k2k],mean)+apply.quarterly(UC[k2k],mean))
SP5.result <- merge(residuals(result.gpuc),predict(result.gpuc),residuals(result.eps),predict(result.eps))

GSPC.predict <- merge(to.monthly(GSPC)[substr(k2k,11,23)],last(spline(seq(1,length(SP5.result[,1]),1),as.vector(SP5.result[,2]),n=length(SP5.result[,1])*3+1)$y,n=length(to.monthly(GSPC)[,1][substr(k2k,11,23)])),last(spline(seq(1,length(SP5.result[,1]),1),as.vector(SP5.result[,4]),n=length(SP5.result[,1])*3+1)$y,n=length(to.monthly(GSPC)[,1][substr(k2k,11,23)])),suffixes=c('','spline','eps'))


plot(merge(GSPC.predict[,4],GSPC.predict[,7],GSPC.predict[,8],GSPC.predict[,4]-GSPC.predict[,7],GSPC.predict[,4]-GSPC.predict[,8]),main="GSPC.predict[,4] vs. GSPC.predict[,7]",grid.ticks.on='months')
tmp.legend <- "Black: actual \nRed: spline\nGreen: eps"
addLegend(legend.loc = "topleft", legend.names = tmp.legend,col=3)
tmp.addTA <- as.xts(rep(2800,length(index(GSPC.predict))),index(GSPC.predict))

addSeries(tmp.addTA,on=1,col=6,lwd=1)

Prepare Data --- getsymbols(), auto.arima(), as.yearmon(), as.yearqtr()

> len_mon <- 48  # # of months to predict
> r <- 1.04 # pesumed GDP growth rate
> i <- seq(2,len_mon/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.
> as.xts(forecast(auto.arima(CS["2012::"]),h=len_mon)$mean[1:len_mon],as.Date(as.yearmon(mondate(index(last(CS)))+seq(1,len_mon,1)),frac=0))[as.Date(as.yearqtr(mondate(index(last(CS)))+seq(3,len_mon,3)),frac=0)]
               [,1]
2018-07-01 225.2251
2018-10-01 227.3579
<skip>
2022-04-01 266.0765
> last(CS)
           SPCS10RSA
2018-06-01  224.8054
> as.xts(forecast(auto.arima(PA),h=len_mon)$mean[1:len_mon],as.Date(as.yearmon(mondate(index(last(PA)))+seq(1,len_mon,1)),frac=0))[as.Date(as.yearqtr(mondate(index(last(PA)))+seq(3,len_mon,3)),frac=0)]
               [,1]
2018-10-01 149684.2
2019-01-01 150292.1
<skip>
2022-07-01 158802.2
>  last(PA)
           PAYEMS
2018-08-01 149279
> last(UC)
           UNDCONTSA
2018-07-01      1122
>
as.xts(forecast(auto.arima(UC),h=len_mon)$mean[1:len_mon],as.Date(as.yearmon(mondate(index(last(UC)))+seq(1,len_mon,1)),frac=0))[as.Date(as.yearqtr(mondate(index(last(UC)))+seq(3,len_mon,3)),frac=0)]
2018-10-01 1130.871
2019-01-01 1139.742
<skip>
2022-07-01 1263.939


len_mon <- 48  # # of months to predict
gdp_g_r <- 1.04 # pesumed GDP growth rate
i <- seq(2,len_mon/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.

getSymbols("PAYEMS",src="FRED")
PA <- PAYEMS
m_PA <- as.xts(forecast(auto.arima(PA),h=len_mon)$mean[1:len_mon],as.Date(as.yearmon(mondate(index(last(PA)))+seq(1,len_mon,1)),frac=0))[as.Date(as.yearqtr(mondate(index(last(PA)))+seq(3,len_mon,3)),frac=0)]
PAq <- apply.quarterly(PA[k2k],mean)
length(PAq)

getSymbols("UNDCONTSA",src="FRED")
UC <- UNDCONTSA
m_UC <- as.xts(forecast(auto.arima(UC),h=len_mon)$mean[1:len_mon],as.Date(as.yearmon(mondate(index(last(UC)))+seq(1,len_mon,1)),frac=0))[as.Date(as.yearqtr(mondate(index(last(UC)))+seq(3,len_mon,3)),frac=0)]
UCq <- apply.quarterly(UC[k2k],mean)
length(UCq)

getSymbols('SPCS10RSA',src='FRED')
CS <- SPCS10RSA
m_CS_2012 <- as.xts(forecast(auto.arima(CS["2012::"]),h=len_mon)$mean[1:len_mon],as.Date(as.yearmon(mondate(index(last(CS)))+seq(1,len_mon,1)),frac=0))[as.Date(as.yearqtr(mondate(index(last(CS)))+seq(3,len_mon,3)),frac=0)]
CSq <- apply.quarterly(CS[k2k],mean)
length(CSq)

getSymbols("GDP",src="FRED")
G <- GDP
# m_GDP <- as.xts(as.vector(last(GDP)) * r**(i/4),d)
m_GDP <- as.xts(as.vector(last(G)) * gdp_g_r**(seq(1,len_mon/3,1)/4),as.Date(as.yearqtr(mondate(index(last(G)))+seq(3,len_mon,3)),frac=0))
kikan <- paste("1992-01-01::",as.Date(as.yearmon((mondate(index(last(G)))+2)),frac=1),sep="")
k2k <- paste("2000-01-01::",as.Date(as.yearmon((mondate(index(last(G)))+2)),frac=1),sep="")

SP5 <- as.xts(read.zoo(read.csv("~/SP5.csv")))

length(CSq)
length(UCq)
length(PAq)
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]))
my_sp5cs(k2k,m_GDP[d[1:9]],m_PA[d[1:9]],m_UC[d[1:9]],m_CS_2012[d[1:9]])
result.eps <- lm(apply.quarterly(SP5[,4][k2k],mean) ~ eps_year_xts[k2k]+apply.quarterly(PA[k2k],mean)+apply.quarterly(CS[k2k],mean)+apply.quarterly(UC[k2k],mean))
result.gpuc <- 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]))
SP5.result <- merge(residuals(result.gpuc),predict(result.gpuc),residuals(result.eps),predict(result.eps))

GSPC.predict <- merge(to.monthly(GSPC)[substr(k2k,11,23)],last(spline(seq(1,length(SP5.result[,1]),1),as.vector(SP5.result[,2]),n=length(SP5.result[,1])*3+1)$y,n=length(to.monthly(GSPC)[,1][substr(k2k,11,23)])),last(spline(seq(1,length(SP5.result[,1]),1),as.vector(SP5.result[,4]),n=length(SP5.result[,1])*3+1)$y,n=length(to.monthly(GSPC)[,1][substr(k2k,11,23)])),suffixes=c('','spline','eps'))

Prepare Data and set parameters.

kg <- "1992-01-01::2018-03-31"
kikan <- "1992-01-01::2018-03-31"
k2k <- "2000-01-01::2018-03-31"   #GPUC model is based data only from 2000 because of cs index availability
# prepare parameters.
l <- 48  # # of months to predict
r <- 1.04 # pesumed GDP growth rate
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.

getSymbols("GDP",src="FRED")
G <- GDP
m_GDP <- as.xts(as.vector(last(GDP)) * r**(i/4),d)

getSymbols("PAYEMS",src="FRED")
PA <- PAYEMS
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]
PAq <- apply.quarterly(PA[k2k],mean)

getSymbols('SPCS10RSA',src='FRED')
CS <- SPCS10RSA
m_CS <-
(as.xts(forecast(auto.arima(CS),h=l)$mean[1:l],as.Date(as.yearmon(seq(mondate(index(last(CS)))+1,by=1,length.out=l))))[(3-month(index(last(CS))) %% 3) + seq(1,l-3,3)])[d]
CSq <- apply.quarterly(CS[k2k],mean)

getSymbols("UNDCONTSA",src="FRED")
# UC <- UNDCONTSA
# comment out here until UNDCONTSA update is resumed in FRED.
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]
UCq <- apply.quarterly(UC[k2k],mean)

How to treat residuals() output as xts object.

> 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]))

above outpus like

2000-03-01   2000-06-01   2000-09-01   2000-12-01   2001-03-01   2001-06-01   2001-09-01
-18.6280839  120.6122215  186.8549495  162.9236931   63.9004346  -15.8557805  -34.5815455
2001-12-01   2002-03-01   2002-06-01   2002-09-01   2002-12-01   2003-03-01   2003-06-01 (....)

and class are like,

> class(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] "xts" "zoo"

however, this is not able to be treated as usual xts object. then do as below.

> resi <- 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]))
> resi[,1]
                   [,1]
2000-03-01  -18.6280839
2000-06-01  120.6122215
2000-09-01  186.8549495
2000-12-01  162.9236931
2001-03-01   63.9004346
(....)

> class(resi[,1])
[1] "xts" "zoo"
>plot(resi[,1])

S&P 500 for next 9 qtrs with and w/o case shiller price index for till 2018Q1

# set the periods of calculation
kg <- "1992-01-01::2018-03-31"
kikan <- "1992-01-01::2018-03-31"
k2k <- "2000-01-01::2018-03-31"   #GPUC model is based data only from 2000 because of cs index availability
# prepare parameters.
l <- 48  # # of months to predict
r <- 1.05 # pesumed GDP growth rate
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.
# for the case the current time is more than 6 months away from the last GDP data, the below is to adjust # of months for GDP forecast.
# get data from FRED and store selfregression data

getSymbols("GDP",src="FRED")
G <- GDP
m_GDP <- as.xts(as.vector(last(GDP)) * r**(i/4),d)
getSymbols("PAYEMS",src="FRED")
PA <- PAYEMS
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]
getSymbols("UNDCONTSA",src="FRED")
# UC <- UNDCONTSA
# comment out here until UNDCONTSA update is resumed in FRED.
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]
getSymbols('SPCS10RSA',src='FRED')
CS <- SPCS10RSA
m_CS <-
(as.xts(forecast(auto.arima(CS),h=l)$mean[1:l],as.Date(as.yearmon(seq(mondate(index(last(CS)))+1,by=1,length.out=l))))[(3-month(index(last(CS))) %% 3) + seq(1,l-3,3)])[d]


# download data from yahoo finace before this.
SP5 <- as.xts(read.zoo(read.csv("~/SP5.csv")))



### caution still in debug
if(floor(MonthsBetween(mondate(last(index(G))),mondate(Sys.Date()))) > 5){ i <- i+1}
### caution ends





# predict next 9 qtrs. by GPU model
my_sp5(kikan,m_GDP[1:9],m_PA_backup[1:9],m_UC[1:9])
# by GPUC model borrow index from PA for case shiller index.
# use "as.xts()", not "xts()"
# this assumes case shiller index of the next qtr is 223 and it will add 2 pts. each qtr.
my_sp5cs(k2k,m_GDP[1:9],m_PA[1:9],m_UC[1:9],as.xts(seq(223,223+2*8,2),index(m_PA[1:9])))
# or use m_CS[1:9]

Draw the graph - Actual, Theory and Residual.

plot(merge(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]))),
+ 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]))))
axis(2,at=c(100,50,0,-50,-100,-150))
axis(4,at=c(100,50,0,-50,-100,-150))

black = actual
red = theory
green = residual

# or draw horizontal line at the designated x-value, please not that
# the index must be synchronized with the main graph, which is equal to
#  seq(as.Date("2000-03-01"),as.Date("2017-12-01"),by='quarters'))
# no to k2k.

plot(merge(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]))),
+            + 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]))),ylim=c(-200,3000))
# addSeries(as.xts(rep(-100,72),index(apply.quarterly(SP5[k2k],mean)[,1])),on=1,ylim=c(-200,2700),lwd=1)
addSeries(as.xts(rep(-100,length(index(apply.quarterly(SP5[k2k],mean)[,1]))),index(apply.quarterly(SP5[k2k],mean)[,1])),on=1,ylim=c(-200,3000),lwd=1)

S&P 500 for next 9 qtrs with and w/o case shiller price index

kg <- "1992-01-01::2017-12-31"
kikan <- "1992-01-01::2017-12-31"
k2k <- "2000-01-01::2017-12-31"   #GPUC model is based data only from 2000 because of cs index availability
# prepare parameters.
l <- 48  # # of months to predict
r <- 1.05 # pesumed GDP growth rate
# get data from FRED
getSymbols("GDP",src="FRED")
getSymbols("PAYEMS",src="FRED")
getSymbols("UNDCONTSA",src="FRED")
getSymbols('SPCS10RSA',src='FRED')
CS <- SPCS10RSA
PA <- PAYEMS
UC <- UNDCONTSA
# download data from yahoo finace before this.
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.
# for the case the current time is more than 6 months away from the last GDP data, the below is to adjust # of months for GDP forecast.
### caution still in debug
if(floor(MonthsBetween(mondate(last(index(G))),mondate(Sys.Date()))) > 5){ i <- i+1}
### caution ends
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]
m_CS <-
(as.xts(forecast(auto.arima(CS),h=l)$mean[1:l],as.Date(as.yearmon(seq(mondate(index(last(CS)))+1,by=1,length.out=l))))[(3-month(index(last(CS))) %% 3) + seq(1,l-3,3)])[d]

# predict next 9 qtrs. by GPU model
my_sp5(kikan,m_GDP[1:9],m_PA_backup[1:9],m_UC[1:9])
# by GPUC model borrow index from PA for case shiller index.
# use "as.xts()", not "xts()"
# this assumes case shiller index of the next qtr is 223 and it will add 2 pts. each qtr.
my_sp5cs(k2k,m_GDP[1:9],m_PA[1:9],m_UC[1:9],as.xts(seq(223,223+2*8,2),index(m_PA[1:9])))
# or use m_CS[1:9]

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