2017年11月2日木曜日

2017Q3 GDP vol1

> library(beepr)
> getSymbols("GDP",src="FRED")
    As of 0.4-0, ‘getSymbols’ uses env=parent.frame() and
 auto.assign=TRUE by default.

 This  behavior  will be  phased out in 0.5-0  when the call  will
 default to use auto.assign=FALSE. getOption("getSymbols.env") and
 getOptions("getSymbols.auto.assign") are now checked for alternate defaults

 This message is shown once per session and may be disabled by setting
 options("getSymbols.warning4.0"=FALSE). See ?getSymbols for more details.
[1] "GDP"
> getSymbols("GDPC96",src="FRED")
[1] "GDPC96"
> getSymbols("PAYEMS",src="FRED")
[1] "PAYEMS"
> getSymbols("UNDCONTSA",src="FRED")
[1] "UNDCONTSA"
> beep(1)
> PA <- PAYEMS
> UC <- UNDCONTSA
> kikan <- "1992-01-01::2017-09-30"
> SP5 <- as.xts(read.zoo(read.csv("~/SP5.csv")))
> G <- GDP
> 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) ))

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

Residuals:
    Min      1Q  Median      3Q     Max
-234.23  -64.92   10.55   75.37  182.71

Coefficients:
                                                                             Estimate Std. Error t value
(Intercept)                                                                 7.671e+03  1.268e+03   6.051
apply.quarterly(PA[kikan], mean)                                           -8.665e-02  1.589e-02  -5.453
G[kikan]                                                                   -8.382e-01  1.593e-01  -5.262
apply.quarterly(PA[kikan], mean):apply.quarterly(UC[kikan], mean)           5.004e-05  7.757e-06   6.451
apply.quarterly(PA[kikan], mean):G[kikan]                                   9.093e-06  1.177e-06   7.724
apply.quarterly(UC[kikan], mean):G[kikan]                                  -7.809e-04  2.113e-04  -3.695
apply.quarterly(PA[kikan], mean):apply.quarterly(UC[kikan], mean):G[kikan]  2.061e-09  1.204e-09   1.711
                                                                           Pr(>|t|)   
(Intercept)                                                                2.77e-08 ***
apply.quarterly(PA[kikan], mean)                                           3.84e-07 ***
G[kikan]                                                                   8.67e-07 ***
apply.quarterly(PA[kikan], mean):apply.quarterly(UC[kikan], mean)          4.49e-09 ***
apply.quarterly(PA[kikan], mean):G[kikan]                                  1.09e-11 ***
apply.quarterly(UC[kikan], mean):G[kikan]                                  0.000365 ***
apply.quarterly(PA[kikan], mean):apply.quarterly(UC[kikan], mean):G[kikan] 0.090318 . 
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 103.5 on 96 degrees of freedom
Multiple R-squared:  0.9594, Adjusted R-squared:  0.9569
F-statistic: 378.5 on 6 and 96 DF,  p-value: < 2.2e-16

> as.xts(forecast(auto.arima(UC),h=10)$mean[1:10],as.Date(as.yearmon(seq(mondate(index(last(UC)))+1,by=1,length.out=10))))[(3-month(index(last(UC))) %% 3) + c(1,4,7)]
               [,1]
2018-01-01 1101.449
2018-04-01 1112.170
2018-07-01 1120.420
> as.xts(forecast(auto.arima(PA),h=10)$mean[1:10],as.Date(as.yearmon(seq(mondate(index(last(PA)))+1,by=1,length.out=10))))[(3-month(index(last(UC))) %% 3) + c(1,4,7)]
               [,1]
2018-01-01 147033.0
2018-04-01 147330.5
2018-07-01 147648.8
> l <- 120  # # of months to predict
> 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.
> r <- 1.05 # pesumed GDP growth rate
> last(PA);last(UC);last(GDP)
           PAYEMS
2017-09-01 146659
           UNDCONTSA
2017-09-01      1082
                GDP
2017-07-01 19495.48
> 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(m_GDP,m_PA,m_UC)
[1] "m_m params! apply.quarter - UC w/ nominal GDP"
               [,1]
2018-01-01 2486.422
2018-04-01 2514.327
2018-07-01 2551.937
2018-10-01 2597.523
2019-01-01 2649.756
2019-04-01 2707.624
<skip>
2026-10-01 5946.552
2027-01-01 6103.463
2027-04-01 6264.051
2027-07-01 6428.380
kg = [1] "1992-01-01::2017-06-30"
> m_PA
               [,1]
2018-01-01 147033.0
2018-04-01 147330.5
2018-07-01 147648.8
2018-10-01 147981.8
2019-01-01 148325.4
2019-04-01 148676.6
<skip>
2026-10-01 159734.3
2027-01-01 160104.5
2027-04-01 160474.7
2027-07-01 160844.9

0 件のコメント: