Calculate day dependent data with less workload - 1 set languge locale

In order to calculate day dependent data with less workload, I like to automate processes. The first step is to judge the first day of each belongs to which day of the week. weekday function provides this function.

Currently R's language setting is set to JA. This may be

> weekdays(as.Date("2016-01-01"),abbreviate = TRUE)
[1] "金"
> Sys.getlocale("LC_MESSAGES")
[1] "ja_JP.UTF-8"

Change LC_MESSAGES to US. But weekdays function still return the value in Japanese.

> Sys.setlocale("LC_MESSAGES",'en_US')
[1] "en_US"
> weekdays(as.Date("2016-01-01"),abbreviate = TRUE)
[1] "金"

Investigate Sys.setlocale and found that there are other classes than LC_MESSAGES. I have no problem with US English.

> Sys.setlocale
function (category = "LC_ALL", locale = "") 
{
    category <- match(category, c("LC_ALL", "LC_COLLATE", "LC_CTYPE", 
        "LC_MONETARY", "LC_NUMERIC", "LC_TIME", "LC_MESSAGES", 
        "LC_PAPER", "LC_MEASUREMENT"))
    if (is.na(category)) 
        stop("invalid 'category' argument")
    .Internal(Sys.setlocale(category, locale))
}
<bytecode: 0x102cfdee0>
<environment: namespace:base>

Set LC_ALL locale, instead of LC_MESSAGES to US English.

> Sys.setlocale("LC_ALL",'en_US')
[1] "en_US/en_US/en_US/C/en_US/en_US"

Then return is in US English.

> weekdays(as.Date("2016-01-01"),abbreviate = TRUE)
[1] "Fri"
>

S&P 500 according to GDP forecast by GDPnow 2016-05-28

Now the target value at the end of 2016 is 2339.

> 0.603*(predict(VAR(v_XCPG_1992_2016,lag=6))$fcst$GDP[,1])-7791
 [1] 2249.308 2266.049 2280.233 2294.785 2306.989 2318.216 2328.660 2339.162 2349.197 2358.947

Set Operator - count working days in 2016

Combined previous two samples will bring the list which are weekday bu not holiday.

> seq(1,365,1)[is.element(seq(1,365,1)%%7,c(6,4,5,0,1))][!is.element(seq(1,365,1)[is.element(seq(1,365,1)%%7,c(6,4,5,0,1))],holiday2016)]


Then you now know how many days you have to go to the office.

> length(seq(1,365,1)[is.element(seq(1,365,1)%%7,c(6,4,5,0,1))][!is.element(seq(1,365,1)[is.element(seq(1,365,1)%%7,c(6,4,5,0,1))],holiday2016)])
[1] 245
> 

Set operator - how to create the list of Monday to Friday seq number for 2016.

Create the list of seq number of weekday for 2016. Here 4 represents Monday. 5 Tuesday, 6 Thursday and 0 for Friday.

> seq(1,365,1)[is.element(seq(1,365,1)%%7,c(6,4,5,0,1))]
  [1]   1   4   5   6   7   8  11  12  13  14  15  18  19  20  21  22  25  26  27  28  29  32  33  34  35  36  39  40  41
 [30]  42  43  46  47  48  49  50  53  54  55  56  57  60  61  62  63  64  67  68  69  70  71  74  75  76  77  78  81  82
 [59]  83  84  85  88  89  90  91  92  95  96  97  98  99 102 103 104 105 106 109 110 111 112 113 116 117 118 119 120 123
 [88] 124 125 126 127 130 131 132 133 134 137 138 139 140 141 144 145 146 147 148 151 152 153 154 155 158 159 160 161 162
[117] 165 166 167 168 169 172 173 174 175 176 179 180 181 182 183 186 187 188 189 190 193 194 195 196 197 200 201 202 203
[146] 204 207 208 209 210 211 214 215 216 217 218 221 222 223 224 225 228 229 230 231 232 235 236 237 238 239 242 243 244
[175] 245 246 249 250 251 252 253 256 257 258 259 260 263 264 265 266 267 270 271 272 273 274 277 278 279 280 281 284 285
[204] 286 287 288 291 292 293 294 295 298 299 300 301 302 305 306 307 308 309 312 313 314 315 316 319 320 321 322 323 326
[233] 327 328 329 330 333 334 335 336 337 340 341 342 343 344 347 348 349 350 351 354 355 356 357 358 361 362 363 364 365
> 

Set Operator - how to extract only working day from all Friday of 2016.

When holiday2016 is the list of seq. number to represent national holidays and seq(1,365,7) is for all Friday in Friday
seq(1,365,7)[!is.element(seq(1,365,7),holiday2016)] will generate the list of sequence numbers for every working Friday in the year.

> holiday2016
 [1]   1   2  11  42  81 120 124 125 126 200 224 263 267 284 308 328 358
> seq(1,365,7)
 [1]   1   8  15  22  29  36  43  50  57  64  71  78  85  92  99 106 113 120 127 134 141 148 155 162 169 176 183 190 197
[30] 204 211 218 225 232 239 246 253 260 267 274 281 288 295 302 309 316 323 330 337 344 351 358 365
> seq(1,365,7)[!is.element(seq(1,365,7),holiday2016)]
 [1]   8  15  22  29  36  43  50  57  64  71  78  85  92  99 106 113 127 134 141 148 155 162 169 176 183 190 197 204 211
[30] 218 225 232 239 246 253 260 274 281 288 295 302 309 316 323 330 337 344 351 365

re-re-forecast GDP

Based upon new update GDPnow forecast at May 28, do re-re-forecast GDP for next 10 months. 2nd qtr. New GDP annual growth rate is 2.8%.

> predict(VAR(v_XCPG_1992_2016,lag.max = 6))$fcst$GDP
          fcst    lower    upper        CI
 [1,] 16650.59 16620.69 16680.49  29.90123
 [2,] 16678.36 16613.26 16743.45  65.09441
 [3,] 16701.88 16596.66 16807.10 105.22027
 [4,] 16726.01 16586.58 16865.44 139.42727
 [5,] 16746.25 16575.42 16917.08 170.83343
 [6,] 16764.87 16565.62 16964.12 199.25079
 [7,] 16782.19 16554.59 17009.78 227.59567
 [8,] 16799.61 16544.35 17054.86 255.25926
 [9,] 16816.25 16533.71 17098.78 282.53532
[10,] 16832.42 16523.43 17141.40 308.98703

Replace a specific column in the data structure.

Prepare updated data in "w". v_XCPG_1992_2016 only has data from 1992. Use "w["1992::"]"  to limit data just from 1992. Data were updated according to GDPnow 2016/5/28.


> tail(v_XCPG_1992_2016$GDP)
2015-11-01 16477.960
2015-12-01 16485.354
2016-01-01 16505.700
2016-02-01 16543.000
2016-03-01 16576.000
2016-04-01 16607.000
> tail(w)
2015-11-01 16482.07
2015-12-01 16493.57
2016-01-01 16505.07
2016-02-01 16540.61
2016-03-01 16580.61
2016-04-01 16620.61
> v_XCPG_1992_2016$GDP <- w["1992::"]
> tail(v_XCPG_1992_2016$GDP)
                GDP
2015-11-01 16482.07
2015-12-01 16493.57
2016-01-01 16505.07
2016-02-01 16540.61
2016-03-01 16580.61
2016-04-01 16620.61

S&P500 forecast next 10 months

With newly updated GDP data and model at  May 25th , S&P500 forecast next 10 months are as below.

Please note that the current price is already more than 100 pts below theoretical value. 2323.292 is expected  at the end of this year.


> 0.603*(predict(VAR(v_XCPG_1992_2016,lag=6))$fcst$GDP[,1])-7791
 [1] 2238.212 2254.216 2268.441 2282.354 2293.338 2303.527 2313.299 2323.292 2332.697 2341.712

GDP forecast by newly updated data.

2016 1st Quarter GDP was revised upward. The growth rate was updated from 0.5% to 0.8%.
The new GDP forecast are as below. The previous one is available at May 23rd entry. #10 months is up 6.1 trillion #8 is also up 6.7. It will also give S&P 500 from another 30 to 40 pts at the end of this year.



> predict(VAR(v_XCPG_1992_2016,lag=6))$fcst$GDP
          fcst    lower    upper        CI
 [1,] 16632.19 16602.29 16662.10  29.90217
 [2,] 16658.73 16593.67 16723.80  65.06516
 [3,] 16682.32 16577.21 16787.44 105.11452
 [4,] 16705.40 16566.14 16844.65 139.25646
 [5,] 16723.61 16552.99 16894.24 170.62578
 [6,] 16740.51 16541.47 16939.54 199.03576
 [7,] 16756.71 16529.33 16984.10 227.38612
 [8,] 16773.29 16518.24 17028.34 255.05082
 [9,] 16788.88 16506.56 17071.21 282.32406
[10,] 16803.83 16495.07 17112.60 308.76632

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

do not confuse with "plot.mlm is not implemented yet"

Please note that lm output class differs according to a given formula. In most cases, it just belongs to "lm". However, it seems that when the target parameter has multiple culumns, result output also belongs to "mlm"

> class(lm(to.monthly(SP500["2009::2016-02-01"]) ~ GDP_XTS["2009::2016-02-01"]))
[1] "mlm" "lm" 
> plot(lm(to.monthly(SP500["2009::2016-02-01"]) ~ GDP_XTS["2009::2016-02-01"]))
 エラー: 'plot.mlm' is not implemented yet
> plot(lm(to.monthly(SP500["2009::2016-02-01"])[,6] ~ GDP_XTS["2009::2016-02-01"]))
> class(lm(to.monthly(SP500["2009::2016-02-01"])[,6] ~ GDP_XTS["2009::2016-02-01"]))
[1] "lm"
> plot(lm(to.monthly(SP500["2009::2016-02-01"])[,6] ~ GDP_XTS["2009::2016-02-01"]))
>

lm(to.quarterly(SP5["2010-04-01::2016-03-01"])[,6] ~ GDPC96["2010-04-01::2016-01-01"] )

A strong correlation exists on those periods, however, the degree of correlation is lower in different terms. It seems that the relationship between GDP and  stock price is more likely to become solid in an upward trend.

> summary(lm(to.quarterly(SP5["2010-04-01::2016-03-01"])[,4] ~ GDPC96["2010-04-01::2016-01-01"] ) )

Call:
lm(formula = to.quarterly(SP5["2010-04-01::2016-03-01"])[, 4] ~  GDPC96["2010-04-01::2016-01-01"])

Residuals:
    Min      1Q  Median      3Q     Max 
-196.63  -85.43   12.53   71.48  171.59 

Coefficients:
                                                      Estimate       Std. Error     t value Pr(>|t|)    
(Intercept)                                       -7.791e+03  6.168e+02  -12.63 1.47e-11 ***
GDPC96["2010-04-01::2016-01-01"]  6.036e-01  3.954e-02   15.27 3.45e-13 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 105.7 on 22 degrees of freedom
Multiple R-squared:  0.9137, Adjusted R-squared:  0.9098 
F-statistic:   233 on 1 and 22 DF,  p-value: 3.453e-13

This model indicates that S&P 500 will up 248 points when GDP annual growth rate is 2.5%.

> last(as.numeric(GDPC96)) * 0.025 * 0.603

[1] 248.6282

Fear for low birth rate and hyper inflation and Gereman balanced budget policy

Fear for low birth rate and hyper inflation and Gereman balanced budget policy. Previously thought that German long-going and unchallenged balanced policy came from hyper inflation after WW1, but now need necessity to see the relationship with low birthrate at those periods. It makes a dramatic contrast with Japan's failure to maintain the policy after "Oil Shock and stagnation" while they didn't have experienced low birthrate in inter-war  period.

calculate X months after the specific date.

"seq" function included in mondate library is enhanced to calculate the distance of months between two dates. This is not possible in the standard "seq" function. Please don't forget to load mondate library beforehand.


> seq(mondate("2009/12/31"), by=1, length.out=8)[8]
mondate: timeunits="months"
[1] 2010/07/31