Alfa Investor In A Good Mood

Nikkei225

2021-05-21T08:39:00.002+09:00

28000-28550 up in the early session, down lately.

Nikkei 225

2021-05-20T08:28:00.004+09:00

27500-28200 expected down in the morning session and up in the afternoon session.

nikkei 225

2021-04-12T13:11:00.004+09:00

今日の予想　３００００円ー２９７８０円

Return Calculus R

2015-06-20T13:37:00.003+09:00

price4<- span=""> read.table('price4.txt',header = T)
price4[1:3,]

##   x5202 x7272 x4927 x4502
## 1    55   682  2460  3595
## 2    56   672  2543  3575
## 3    56   668  2544  3570

summary(price4)

##      x5202           x7272            x4927          x4502     
##  Min.   : 55.0   Min.   : 668.0   Min.   :2363   Min.   :3500  
##  1st Qu.: 73.0   1st Qu.: 780.0   1st Qu.:2467   1st Qu.:3710  
##  Median : 96.0   Median : 948.0   Median :2525   Median :3865  
##  Mean   : 92.4   Mean   : 968.8   Mean   :2633   Mean   :4199  
##  3rd Qu.:110.0   3rd Qu.:1098.0   3rd Qu.:2832   3rd Qu.:4700  
##  Max.   :124.0   Max.   :1334.0   Max.   :3085   Max.   :5280

NipponSheetGlassCompany_TSE5202<- span=""> price4$x5202
PolaOrbisHoldings_TSE4927<- span=""> price4$x4927
YamahaMotor_TSE7272<- span=""> price4$x7272
TakedaPharmaceuticalCompany_TSE4502<- span=""> price4$x4502
par(mfcol=c(2,2))
plot(NipponSheetGlassCompany_TSE5202,type="l"
     ,main="NipponSheetGlassCompany_TSE5202")
plot(PolaOrbisHoldings_TSE4927,type="l"
     ,main="PolaOrbisHoldings_TSE4927")
plot(YamahaMotor_TSE7272,type="l"
     ,main="YamahaMotor_TSE7272")
plot(TakedaPharmaceuticalCompany_TSE4502,type="l"
     ,main="TakedaPharmaceuticalCompany_TSE4502")

時系列の始点を基準として、価格収益率を求める

par(mfcol=c(2,2))
plot(NipponSheetGlassCompany_TSE5202/NipponSheetGlassCompany_TSE5202[1]
     ,type="l",main="NipponSheetGlassCompany_TSE5202")
plot(PolaOrbisHoldings_TSE4927/PolaOrbisHoldings_TSE4927[1]
     ,type="l",main="PolaOrbisHoldings_TSE4927")
plot(YamahaMotor_TSE7272/YamahaMotor_TSE7272[1]
     ,type="l",main="YamahaMotor_TSE7272")
plot(TakedaPharmaceuticalCompany_TSE4502/TakedaPharmaceuticalCompany_TSE4502[1]
     ,type="l",main="TakedaPharmaceuticalCompany_TSE4502")

判明した点板硝子は最高で2.2倍に株価が成長した

次に対数差収益を求める##

par(mfcol=c(2,2))
plot(diff(log(NipponSheetGlassCompany_TSE5202)*100)
     ,type="l",main="NipponSheetGlassCompany_TSE5202"
     ,xlab="time",ylab="NipponSheetGlassCompany_TSE5202")
plot(diff(log(PolaOrbisHoldings_TSE4927)*100)
     ,type="l",main="PolaOrbisHoldings_TSE4927"
     ,xlab="time",ylab="NipponSheetGlassCompany_TSE5202")
plot(diff(log(YamahaMotor_TSE7272)*100)
     ,type="l",main="YamahaMotor_TSE7272"
     ,xlab="time",ylab="NipponSheetGlassCompany_TSE5202")
plot(diff(log(TakedaPharmaceuticalCompany_TSE4502)*100)
     ,type="l",main="TakedaPharmaceuticalCompany_TSE4502"
     ,xlab="time",ylab="NipponSheetGlassCompany_TSE5202")

結論板硝子とYamahaMotor_TSE7272のボラは高い。一方、武田薬品のそれは低い

par(mfcol=c(2,2))
plot(diff(log(NipponSheetGlassCompany_TSE5202)*100)
     ,type="l",main="NipponSheetGlassCompany_TSE5202"
     ,xlab="time",ylab="NipponSheetGlassCompany_TSE5202"
     ,ylim=range(diff(log(NipponSheetGlassCompany_TSE5202))*100))
plot(diff(log(PolaOrbisHoldings_TSE4927)*100)
     ,type="l",main="PolaOrbisHoldings_TSE4927"
     ,xlab="time",ylab="PolaOrbisHoldings_TSE4927"
     ,ylim=range(diff(log(PolaOrbisHoldings_TSE4927))*100))
plot(diff(log(YamahaMotor_TSE7272)*100)
     ,type="l",main="YamahaMotor_TSE7272"
     ,xlab="time",ylab="YamahaMotor_TSE7272"
     ,ylim=range(diff(log(YamahaMotor_TSE7272))*100))
plot(diff(log(TakedaPharmaceuticalCompany_TSE4502)*100)
     ,type="l",main="TakedaPharmaceuticalCompany_TSE4502"
     ,xlab="time",ylab="TakedaPharmaceuticalCompany_TSE4502"
     ,ylim=range(diff(log(TakedaPharmaceuticalCompany_TSE4502))*100))

################################## ########時系列データ分布と要約########### ##4銘柄の平均収益率

mean(TakedaPharmaceuticalCompany_TSE4502)

## [1] 4199.38

mean(PolaOrbisHoldings_TSE4927)

## [1] 2633.157

mean(NipponSheetGlassCompany_TSE5202)

## [1] 92.39669

mean(YamahaMotor_TSE7272)

## [1] 968.8264

4銘柄のボラティリティ

sd(NipponSheetGlassCompany_TSE5202)

## [1] 20.12854

sd(PolaOrbisHoldings_TSE4927)

## [1] 209.8629

sd(YamahaMotor_TSE7272)

## [1] 187.3384

sd(TakedaPharmaceuticalCompany_TSE4502)

## [1] 571.4873

μ、σ、histgram

par(mfcol=c(2,2))
hist(NipponSheetGlassCompany_TSE5202,main="NipponSheetGlassCompany_TSE5202")
hist(PolaOrbisHoldings_TSE4927,main="PolaOrbisHoldings_TSE4927")
hist(YamahaMotor_TSE7272,main="YamahaMotor_TSE7272")
hist(TakedaPharmaceuticalCompany_TSE4502,main="TakedaPharmaceuticalCompany_TSE4502")

散布図行列

cor(YamahaMotor_TSE7272,TakedaPharmaceuticalCompany_TSE4502)

## [1] 0.9539355

plot(price4)

cor(price4)##相関係数行列

##           x5202     x7272     x4927     x4502
## x5202 1.0000000 0.8851551 0.7087943 0.8293775
## x7272 0.8851551 1.0000000 0.8705264 0.9539355
## x4927 0.7087943 0.8705264 1.0000000 0.9455436
## x4502 0.8293775 0.9539355 0.9455436 1.0000000

Shapiro-Wilk検定

shapiro.test(PolaOrbisHoldings_TSE4927)

## 
##  Shapiro-Wilk normality test
## 
## data:  PolaOrbisHoldings_TSE4927
## W = 0.85788, p-value = 2.071e-09

shapiro.test(TakedaPharmaceuticalCompany_TSE4502)

## 
##  Shapiro-Wilk normality test
## 
## data:  TakedaPharmaceuticalCompany_TSE4502
## W = 0.86407, p-value = 3.767e-09

shapiro.test(NipponSheetGlassCompany_TSE5202)

## 
##  Shapiro-Wilk normality test
## 
## data:  NipponSheetGlassCompany_TSE5202
## W = 0.91089, p-value = 6.758e-07

shapiro.test(YamahaMotor_TSE7272)

## 
##  Shapiro-Wilk normality test
## 
## data:  YamahaMotor_TSE7272
## W = 0.94524, p-value = 9.178e-05

　連の検定(runs test)

library(tseries)
yamaha=as.factor(YamahaMotor_TSE7272<mean(YamahaMotor_TSE7272))
runs.test(yamaha)

## 
##  Runs Test
## 
## data:  yamaha
## Standard Normal = -10.863, p-value < 2.2e-16
## alternative hypothesis: two.sided

Explore the time dependency for four stocks

par(mfcol=c(2,1))
plot(NipponSheetGlassCompany_TSE5202,type='l',main="NipponSheetGlassCompany_TSE5202")
plot(NipponSheetGlassCompany_TSE5202/NipponSheetGlassCompany_TSE5202[1]
     ,type="l",main="NipponSheetGlassCompany_TSE5202")
abline(h=0,lty=3)

par(mfcol=c(2,1))
plot(NipponSheetGlassCompany_TSE5202,type = 'l',main="PolaOrbisHoldings_TSE4927")
plot(PolaOrbisHoldings_TSE4927/PolaOrbisHoldings_TSE4927[1]
     ,type="l",main="PolaOrbisHoldings_TSE4927")
abline(h=0,lty=3)

par(mfcol=c(2,1))
plot(YamahaMotor_TSE7272,type = 'l',main="YamahaMotor_TSE7272")
plot(YamahaMotor_TSE7272/YamahaMotor_TSE7272[1]
     ,type="l",main="YamahaMotor_TSE7272")
abline(h=0,lty=3)

par(mfcol=c(2,1))
plot(TakedaPharmaceuticalCompany_TSE4502,type = 'l',main="TakedaPharmaceuticalCompany_TSE4502")
plot(TakedaPharmaceuticalCompany_TSE4502/TakedaPharmaceuticalCompany_TSE4502[1]
     ,type="l",main="TakedaPharmaceuticalCompany_TSE4502")
abline(h=0,lty=3)

analyze four historical time series

2015-06-19T21:21:00.002+09:00

In [102]:

import pandas as pd
import pandas.io.data as web
import numpy as np
import numpy.random as npr
import statsmodels.api as sm
import scipy.stats as scs
import matplotlib.pyplot as plt
%matplotlib inline



wrapper function for the describe function from

the scipy.stats sublibrary




In [104]:

def print_statistics(array):
    '''Prints selected statistics.
    Parameters
    array: ndarray
    object to generate statistics on
    '''

    sta=scs.describe(array)
    print ' %14s %15s ' %(' statistics ' ,' value' )
    print 30* '-'
    print ' %14s %15f ' %('size  ' ,sta[0] )
    print ' %14s %15f ' %('min  ' , sta[1][0])
    print ' %14s %15f ' %(' max ' ,sta[1][1] )
    print ' %14s %15f ' %(' mean ' ,sta[2] )
    print ' %14s %15f ' %('  std' ,np.sqrt(sta[3]) )
    print ' %14s %15f ' %(' skew ' ,sta[4] )
    print ' %14s %15f ' %('  kurtosis' , sta[5])




The function normality_tests combines three different statistical tests: 
Skewness test ( skewtest ) This tests whether the skew of the sample data 
is “normal” (i.e., has a value close enough to zero). Kurtosis test ( kurtosistest )

Similarly, this tests whether the kurtosis of the sample data is “normal”

(again, close enough to zero). Normality test ( normaltest ) This combines

the other two test approaches to test for normality.




In [105]:

def normality_tests(arr):
    '''Tests for normality distribution of given data set.
Parameters
array: ndarray
object to generate statistics on
    '''
    print ' Skew of data set %14.3f' %scs.skew(arr)
    print ' Skew test p-value  %14.3f' %scs.skewtest(arr)[1]
    print '  kurt of data set %14.3f' %scs.kurtosis(arr)
    print '  kurt test p-value %14.3f' %scs.kurtosistest(arr)[1]
    print '  Norm test p-values %14.3f' %scs.normaltest(arr)[1]
   





the German DAX index :^GDAXI the American S&P 500 index:

^GSPC YHOO: Yahoo MSFT: Micro Soft




In [106]:

symbols=['^GDAXI','^GSPC','YHOO','MSFT']


In [107]:

symbols


Out[107]:

['^GDAXI', '^GSPC', 'YHOO', 'MSFT']

In [108]:

data=pd.DataFrame()
for sym in symbols:
    data[sym]=web.DataReader(sym,data_source='yahoo',start='1/1/2006')['Adj Close']
data=data.dropna()


In [109]:

data.info()




DatetimeIndex: 2352 entries, 2006-01-03 00:00:00 to 2015-06-15 00:00:00
Data columns (total 4 columns):
^GDAXI    2352 non-null float64
^GSPC     2352 non-null float64
YHOO      2352 non-null float64
MSFT      2352 non-null float64
dtypes: float64(4)



Index start at 100

In [110]:

(data/data.ix[0]*100).plot(figsize=(8,6))


Out[110]:





In [111]:

log_returns=np.log(data/data.shift(1))


In [112]:

log_returns.head()


Out[112]:

^GDAXI^GSPCYHOOMSFT
Date
2006-01-03NaNNaNNaNNaN
2006-01-040.0114600.0036660.0014660.004832
2006-01-05-0.0012840.0000160.0135760.000741
2006-01-060.0035810.0093560.039656-0.002968
2006-01-090.0001430.0036500.004848-0.001860


In [113]:

log_returns.hist(bins=50,figsize=(9,6))


Out[113]:

array([[,
        ],
       [,
        ]], dtype=object)



In [114]:

for sym in symbols:
    print ' \n Results for symbols %s' %sym
    print 30*'-'
    log_data=np.array(log_returns[sym].dropna())
    print_statistics(log_data)



 
 Results for symbols ^GDAXI
------------------------------
    statistics            value 
------------------------------
         size       2351.000000 
          min         -0.077391 
           max         0.107975 
          mean         0.000297 
            std        0.014518 
          skew         0.015973 
       kurtosis        6.204635 
 
 Results for symbols ^GSPC
------------------------------
    statistics            value 
------------------------------
         size       2351.000000 
          min         -0.094695 
           max         0.109572 
          mean         0.000211 
            std        0.013282 
          skew        -0.318551 
       kurtosis       10.343515 
 
 Results for symbols YHOO
------------------------------
    statistics            value 
------------------------------
         size       2351.000000 
          min         -0.246364 
           max         0.391817 
          mean        -0.000005 
            std        0.025672 
          skew         0.544870 
       kurtosis       32.247736 
 
 Results for symbols MSFT
------------------------------
    statistics            value 
------------------------------
         size       2351.000000 
          min         -0.124578 
           max         0.170626 
          mean         0.000316 
            std        0.017844 
          skew         0.047480 
       kurtosis       10.277228 


In [115]:

sm.qqplot(log_returns['^GSPC'].dropna(),line='s')
plt.grid(True)
plt.xlabel('theoretical quantiles')
plt.ylabel('sample quantiles')


Out[115]:





In [116]:

sm.qqplot(log_returns['MSFT'].dropna(),line='s')
plt.grid(True)
plt.xlabel('theoretical quantiles')
plt.ylabel('sample quantiles')


Out[116]:





In [117]:

sm.qqplot(log_returns['^GDAXI'].dropna(),line='s')
plt.grid(True)
plt.xlabel('theoretical quantiles')
plt.ylabel('sample quantiles')


Out[117]:





In [118]:

sm.qqplot(log_returns['YHOO'].dropna(),line='s')
plt.grid(True)
plt.xlabel('theoretical quantiles')
plt.ylabel('sample quantiles')


Out[118]:





In [119]:

for sym in symbols:
    print ' \nResults for symbol %s' %sym
    print 32*'-'
    log_data=np.array(log_returns[sym].dropna())
    normality_tests(log_data)



 
Results for symbol ^GDAXI
--------------------------------
 Skew of data set          0.016
 Skew test p-value           0.751
  kurt of data set          6.205
  kurt test p-value          0.000
  Norm test p-values          0.000
 
Results for symbol ^GSPC
--------------------------------
 Skew of data set         -0.319
 Skew test p-value           0.000
  kurt of data set         10.344
  kurt test p-value          0.000
  Norm test p-values          0.000
 
Results for symbol YHOO
--------------------------------
 Skew of data set          0.545
 Skew test p-value           0.000
  kurt of data set         32.248
  kurt test p-value          0.000
  Norm test p-values          0.000
 
Results for symbol MSFT
--------------------------------
 Skew of data set          0.047
 Skew test p-value           0.346
  kurt of data set         10.277
  kurt test p-value          0.000
  Norm test p-values          0.000

Installed pymc library Python

2015-06-16T23:04:00.004+09:00

pymc 2.3.4 --Python Library(Markov Chain Monte Carlo sampling toolkit)

Go to the download site here .

I am using Ubuntu 14.04
shelling:
mkdir pymc 2.3.4
cd pymc 2.3.4
sudo python setup.py install

Time dependency

2015-06-06T16:31:00.000+09:00

price4<- span=""> read.table('price4.txt',header = T)

Historical price

NipponSheetGlassCompany_TSE5202<- span=""> price4$x5202
PolaOrbisHoldings_TSE4927<- span=""> price4$x4927
YamahaMotor_TSE7272<- span=""> price4$x7272
TakedaPharmaceuticalCompany_TSE4502<- span=""> price4$x4502

Return on price start at time 0

Return_TSE5202<- span=""> NipponSheetGlassCompany_TSE5202/NipponSheetGlassCompany_TSE5202[1]
Return_TSE4927<- span=""> PolaOrbisHoldings_TSE4927/PolaOrbisHoldings_TSE4927[1]
Return_TSE7272<- span=""> YamahaMotor_TSE7272/YamahaMotor_TSE7272[1]
Return_TSE4502<- span=""> TakedaPharmaceuticalCompany_TSE4502/TakedaPharmaceuticalCompany_TSE4502[1]

par(mfcol=c(2,2))
plot(Return_TSE5202[1:120],Return_TSE5202[2:121])
plot(Return_TSE4927[1:120],Return_TSE4927[2:121])
plot(Return_TSE7272[1:120],Return_TSE7272[2:121])
plot(Return_TSE4502[1:120],Return_TSE4502[2:121])

cor(Return_TSE5202[1:120],Return_TSE5202[2:121])

## [1] 0.9853242

cor(Return_TSE4927[1:120],Return_TSE4927[2:121])

## [1] 0.9800802

cor(Return_TSE7272[1:120],Return_TSE7272[2:121])

## [1] 0.9893316

cor(Return_TSE4502[1:120],Return_TSE4502[2:121])

## [1] 0.9960337

acf(Return_TSE5202,plot=F)

## 
## Autocorrelations of series 'Return_TSE5202', by lag
## 
##     0     1     2     3     4     5     6     7     8     9    10    11 
## 1.000 0.969 0.937 0.908 0.885 0.866 0.846 0.822 0.800 0.774 0.750 0.727 
##    12    13    14    15    16    17    18    19    20 
## 0.705 0.685 0.665 0.646 0.629 0.607 0.582 0.554 0.524

acf(Return_TSE5202)

acf(Return_TSE4927,plot=F)

## 
## Autocorrelations of series 'Return_TSE4927', by lag
## 
##     0     1     2     3     4     5     6     7     8     9    10    11 
## 1.000 0.965 0.943 0.913 0.885 0.860 0.835 0.809 0.788 0.767 0.743 0.710 
##    12    13    14    15    16    17    18    19    20 
## 0.684 0.656 0.629 0.602 0.579 0.544 0.526 0.504 0.483

acf(Return_TSE4927)

acf(Return_TSE7272,plot=F)

## 
## Autocorrelations of series 'Return_TSE7272', by lag
## 
##     0     1     2     3     4     5     6     7     8     9    10    11 
## 1.000 0.967 0.937 0.902 0.871 0.840 0.808 0.777 0.744 0.716 0.681 0.650 
##    12    13    14    15    16    17    18    19    20 
## 0.621 0.599 0.577 0.558 0.541 0.524 0.510 0.497 0.483

acf(Return_TSE7272)

acf(Return_TSE4502,plot=F)

## 
## Autocorrelations of series 'Return_TSE4502', by lag
## 
##     0     1     2     3     4     5     6     7     8     9    10    11 
## 1.000 0.982 0.962 0.941 0.915 0.890 0.865 0.839 0.813 0.788 0.759 0.732 
##    12    13    14    15    16    17    18    19    20 
## 0.707 0.680 0.654 0.628 0.604 0.581 0.558 0.535 0.515

acf(Return_TSE4502)

acf(Return_TSE5202,plot=F,type="partial") #or type="p"

## 
## Partial autocorrelations of series 'Return_TSE5202', by lag
## 
##      1      2      3      4      5      6      7      8      9     10 
##  0.969 -0.033  0.033  0.095  0.028  0.002 -0.074  0.039 -0.086  0.014 
##     11     12     13     14     15     16     17     18     19     20 
## -0.010 -0.006  0.032 -0.009  0.014  0.018 -0.079 -0.048 -0.079 -0.061

acf(Return_TSE5202,plot=T,type="partial") #or type="p"

acf(Return_TSE4927,plot=F,type="partial") #or type="p"

## 
## Partial autocorrelations of series 'Return_TSE4927', by lag
## 
##      1      2      3      4      5      6      7      8      9     10 
##  0.965  0.159 -0.079 -0.032  0.037  0.016 -0.034  0.038  0.019 -0.045 
##     11     12     13     14     15     16     17     18     19     20 
## -0.180  0.042  0.028 -0.028 -0.026  0.047 -0.166  0.141  0.028  0.002

acf(Return_TSE4927,plot=T,type="partial") #or type="p"

acf(Return_TSE7272,plot=F,type="partial") #or type="p"

## 
## Partial autocorrelations of series 'Return_TSE7272', by lag
## 
##      1      2      3      4      5      6      7      8      9     10 
##  0.967  0.024 -0.082  0.026 -0.005 -0.038  0.003 -0.044  0.055 -0.122 
##     11     12     13     14     15     16     17     18     19     20 
##  0.022  0.042  0.072 -0.006  0.037  0.002  0.020 -0.006  0.041 -0.039

acf(Return_TSE7272,plot=T,type="partial") #or type="p"

acf(Return_TSE4502,plot=F,type="partial") #or type="p"

## 
## Partial autocorrelations of series 'Return_TSE4502', by lag
## 
##      1      2      3      4      5      6      7      8      9     10 
##  0.982 -0.092 -0.028 -0.116  0.004 -0.005 -0.058  0.031 -0.007 -0.134 
##     11     12     13     14     15     16     17     18     19     20 
##  0.070  0.009 -0.026 -0.032 -0.010  0.083 -0.039 -0.031  0.031  0.013

acf(Return_TSE4502,plot=T,type="partial") #or type="p"

Ljung=Box test for X

Box.test(Return_TSE5202,type="L")

## 
##  Box-Ljung test
## 
## data:  Return_TSE5202
## X-squared = 116.45, df = 1, p-value < 2.2e-16

Box.test(Return_TSE4927,type="L")

## 
##  Box-Ljung test
## 
## data:  Return_TSE4927
## X-squared = 115.57, df = 1, p-value < 2.2e-16

Box.test(Return_TSE7272,type="L")

## 
##  Box-Ljung test
## 
## data:  Return_TSE7272
## X-squared = 116.02, df = 1, p-value < 2.2e-16

Box.test(Return_TSE4502,type="L")

## 
##  Box-Ljung test
## 
## data:  Return_TSE4502
## X-squared = 119.7, df = 1, p-value < 2.2e-16

mean(Return_TSE5202)

## [1] 1.67994

mean(Return_TSE4927)

## [1] 1.070389

mean(Return_TSE7272)

## [1] 1.420567

mean(Return_TSE4502)

## [1] 1.168117

acf(Return_TSE5202,plot=F,type="covariance")

## 
## Autocovariances of series 'Return_TSE5202', by lag
## 
##      0      1      2      3      4      5      6      7      8      9 
## 0.1328 0.1287 0.1244 0.1206 0.1176 0.1150 0.1124 0.1092 0.1063 0.1028 
##     10     11     12     13     14     15     16     17     18     19 
## 0.0997 0.0966 0.0936 0.0910 0.0884 0.0858 0.0835 0.0806 0.0773 0.0736 
##     20 
## 0.0697

acf(Return_TSE4927,plot=F,type="covariance")

## 
## Autocovariances of series 'Return_TSE4927', by lag
## 
##       0       1       2       3       4       5       6       7       8 
## 0.00722 0.00697 0.00680 0.00659 0.00639 0.00620 0.00603 0.00584 0.00569 
##       9      10      11      12      13      14      15      16      17 
## 0.00553 0.00536 0.00513 0.00493 0.00473 0.00454 0.00435 0.00418 0.00393 
##      18      19      20 
## 0.00380 0.00363 0.00349

acf(Return_TSE7272,plot=F,type="covariance")

## 
## Autocovariances of series 'Return_TSE7272', by lag
## 
##      0      1      2      3      4      5      6      7      8      9 
## 0.0748 0.0724 0.0701 0.0675 0.0652 0.0628 0.0604 0.0581 0.0557 0.0536 
##     10     11     12     13     14     15     16     17     18     19 
## 0.0510 0.0486 0.0465 0.0448 0.0432 0.0418 0.0405 0.0392 0.0381 0.0372 
##     20 
## 0.0361

acf(Return_TSE4502,plot=F,type="covariance")

## 
## Autocovariances of series 'Return_TSE4502', by lag
## 
##      0      1      2      3      4      5      6      7      8      9 
## 0.0251 0.0246 0.0241 0.0236 0.0229 0.0223 0.0217 0.0210 0.0204 0.0198 
##     10     11     12     13     14     15     16     17     18     19 
## 0.0190 0.0184 0.0177 0.0170 0.0164 0.0157 0.0151 0.0146 0.0140 0.0134 
##     20 
## 0.0129

ar(Return_TSE5202,aic=F,order.max=1)

## 
## Call:
## ar(x = Return_TSE5202, aic = F, order.max = 1)
## 
## Coefficients:
##     1  
## 0.969  
## 
## Order selected 1  sigma^2 estimated as  0.00825

(1-0.631)*mean(Return_TSE5202) #μを求める

## [1] 0.6198978

ar(Return_TSE5202,aic=T,oreder.max=1)

## 
## Call:
## ar(x = Return_TSE5202, aic = T, oreder.max = 1)
## 
## Coefficients:
##     1  
## 0.969  
## 
## Order selected 1  sigma^2 estimated as  0.00825

ar(Return_TSE4927,aic=F,order.max=1)

## 
## Call:
## ar(x = Return_TSE4927, aic = F, order.max = 1)
## 
## Coefficients:
##      1  
## 0.9653  
## 
## Order selected 1  sigma^2 estimated as  0.0005005

(1-0.631)*mean(Return_TSE4927) #μを求める

## [1] 0.3949736

ar(Return_TSE4927,aic=T,oreder.max=1)

## 
## Call:
## ar(x = Return_TSE4927, aic = T, oreder.max = 1)
## 
## Coefficients:
##      1       2  
## 0.8118  0.1590  
## 
## Order selected 2  sigma^2 estimated as  0.000492

ar(Return_TSE7272,aic=F,order.max=1)

## 
## Call:
## ar(x = Return_TSE7272, aic = F, order.max = 1)
## 
## Coefficients:
##      1  
## 0.9672  
## 
## Order selected 1  sigma^2 estimated as  0.004912

(1-0.631)*mean(Return_TSE7272) #μを求める

## [1] 0.5241891

ar(Return_TSE7272,aic=T,oreder.max=1)

## 
## Call:
## ar(x = Return_TSE7272, aic = T, oreder.max = 1)
## 
## Coefficients:
##      1  
## 0.9672  
## 
## Order selected 1  sigma^2 estimated as  0.004912

ar(Return_TSE4502,aic=F,order.max=1)

## 
## Call:
## ar(x = Return_TSE4502, aic = F, order.max = 1)
## 
## Coefficients:
##      1  
## 0.9824  
## 
## Order selected 1  sigma^2 estimated as  0.0008892

(1-0.631)*mean(Return_TSE4502) #μを求める

## [1] 0.4310351

ar(Return_TSE4502,aic=T,oreder.max=1)

## 
## Call:
## ar(x = Return_TSE4502, aic = T, oreder.max = 1)
## 
## Coefficients:
##      1  
## 0.9824  
## 
## Order selected 1  sigma^2 estimated as  0.0008892

Dickey-Fuller test

library(tseries)
adf.test(NipponSheetGlassCompany_TSE5202)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  NipponSheetGlassCompany_TSE5202
## Dickey-Fuller = -1.0795, Lag order = 4, p-value = 0.9223
## alternative hypothesis: stationary

adf.test(Return_TSE5202)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  Return_TSE5202
## Dickey-Fuller = -1.0795, Lag order = 4, p-value = 0.9223
## alternative hypothesis: stationary

adf.test(PolaOrbisHoldings_TSE4927)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  PolaOrbisHoldings_TSE4927
## Dickey-Fuller = -2.3731, Lag order = 4, p-value = 0.4216
## alternative hypothesis: stationary

adf.test(Return_TSE4927)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  Return_TSE4927
## Dickey-Fuller = -2.3731, Lag order = 4, p-value = 0.4216
## alternative hypothesis: stationary

adf.test(YamahaMotor_TSE7272)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  YamahaMotor_TSE7272
## Dickey-Fuller = -2.9204, Lag order = 4, p-value = 0.1943
## alternative hypothesis: stationary

adf.test(Return_TSE7272)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  Return_TSE7272
## Dickey-Fuller = -2.9204, Lag order = 4, p-value = 0.1943
## alternative hypothesis: stationary

adf.test(TakedaPharmaceuticalCompany_TSE4502)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  TakedaPharmaceuticalCompany_TSE4502
## Dickey-Fuller = -1.9524, Lag order = 4, p-value = 0.5964
## alternative hypothesis: stationary

adf.test(Return_TSE4502)

## 
##  Augmented Dickey-Fuller Test
## 
## data:  Return_TSE4502
## Dickey-Fuller = -1.9524, Lag order = 4, p-value = 0.5964
## alternative hypothesis: stationary

Wine Prediction

2014-12-24T05:05:00.000+09:00

wine<-read .scv="" comment----="">('wine.csv')
str(wine)

## 'data.frame':    25 obs. of  7 variables:
##  $ Year       : int  1952 1953 1955 1957 1958 1959 1960 1961 1962 1963 ...
##  $ Price      : num  7.5 8.04 7.69 6.98 6.78 ...
##  $ WinterRain : int  600 690 502 420 582 485 763 830 697 608 ...
##  $ AGST       : num  17.1 16.7 17.1 16.1 16.4 ...
##  $ HarvestRain: int  160 80 130 110 187 187 290 38 52 155 ...
##  $ Age        : int  31 30 28 26 25 24 23 22 21 20 ...
##  $ FrancePop  : num  43184 43495 44218 45152 45654 ...

summary(wine)

##       Year          Price         WinterRain         AGST      
##  Min.   :1952   Min.   :6.205   Min.   :376.0   Min.   :14.98  
##  1st Qu.:1960   1st Qu.:6.519   1st Qu.:536.0   1st Qu.:16.20  
##  Median :1966   Median :7.121   Median :600.0   Median :16.53  
##  Mean   :1966   Mean   :7.067   Mean   :605.3   Mean   :16.51  
##  3rd Qu.:1972   3rd Qu.:7.495   3rd Qu.:697.0   3rd Qu.:17.07  
##  Max.   :1978   Max.   :8.494   Max.   :830.0   Max.   :17.65  
##   HarvestRain         Age         FrancePop    
##  Min.   : 38.0   Min.   : 5.0   Min.   :43184  
##  1st Qu.: 89.0   1st Qu.:11.0   1st Qu.:46584  
##  Median :130.0   Median :17.0   Median :50255  
##  Mean   :148.6   Mean   :17.2   Mean   :49694  
##  3rd Qu.:187.0   3rd Qu.:23.0   3rd Qu.:52894  
##  Max.   :292.0   Max.   :31.0   Max.   :54602

model1<- class="identifier" color="#000000" font="">lm(Price~AGST,data=wine)
summary(model1)

## 
## Call:
## lm(formula = Price ~ AGST, data = wine)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.78450 -0.23882 -0.03727  0.38992  0.90318 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  -3.4178     2.4935  -1.371 0.183710    
## AGST          0.6351     0.1509   4.208 0.000335 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.4993 on 23 degrees of freedom
## Multiple R-squared:  0.435,  Adjusted R-squared:  0.4105 
## F-statistic: 17.71 on 1 and 23 DF,  p-value: 0.000335

model1$residuals

##           1           2           3           4           5           6 
##  0.04204258  0.82983774  0.21169394  0.15609432 -0.23119140  0.38991701 
##           7           8           9          10          11          12 
## -0.48959140  0.90318115  0.45372410  0.14887461 -0.23882157 -0.08974238 
##          13          14          15          16          17          18 
##  0.66185660 -0.05211511 -0.62726647 -0.74714947  0.42113502 -0.03727441 
##          19          20          21          22          23          24 
##  0.10685278 -0.78450270 -0.64017590 -0.05508720 -0.67055321 -0.22040381 
##          25 
##  0.55866518

SSE<- class="identifier" color="#000000" font="">sum(model1$residuals^2)
SSE

## [1] 5.734875

model2<- class="identifier" color="#000000" font="">lm(Price~AGST+HarvestRain,data=wine)
summary(model2)

## 
## Call:
## lm(formula = Price ~ AGST + HarvestRain, data = wine)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.88321 -0.19600  0.06178  0.15379  0.59722 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -2.20265    1.85443  -1.188 0.247585    
## AGST         0.60262    0.11128   5.415 1.94e-05 ***
## HarvestRain -0.00457    0.00101  -4.525 0.000167 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3674 on 22 degrees of freedom
## Multiple R-squared:  0.7074, Adjusted R-squared:  0.6808 
## F-statistic: 26.59 on 2 and 22 DF,  p-value: 1.347e-06

SSE<- class="identifier" color="#000000" font="">sum(model2$residuals^2)
SSE

## [1] 2.970373

model3<- class="identifier" color="#000000" font="">lm(Price~AGST+HarvestRain+WinterRain+Age+FrancePop,data=wine)
model3

## 
## Call:
## lm(formula = Price ~ AGST + HarvestRain + WinterRain + Age + 
##     FrancePop, data = wine)
## 
## Coefficients:
## (Intercept)         AGST  HarvestRain   WinterRain          Age  
##  -4.504e-01    6.012e-01   -3.958e-03    1.043e-03    5.847e-04  
##   FrancePop  
##  -4.953e-05

SSE<- class="identifier" color="#000000" font="">sum(model3$residuals^2)
SSE

## [1] 1.732113

summary(model3)

## 
## Call:
## lm(formula = Price ~ AGST + HarvestRain + WinterRain + Age + 
##     FrancePop, data = wine)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.48179 -0.24662 -0.00726  0.22012  0.51987 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -4.504e-01  1.019e+01  -0.044 0.965202    
## AGST         6.012e-01  1.030e-01   5.836 1.27e-05 ***
## HarvestRain -3.958e-03  8.751e-04  -4.523 0.000233 ***
## WinterRain   1.043e-03  5.310e-04   1.963 0.064416 .  
## Age          5.847e-04  7.900e-02   0.007 0.994172    
## FrancePop   -4.953e-05  1.667e-04  -0.297 0.769578    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.3019 on 19 degrees of freedom
## Multiple R-squared:  0.8294, Adjusted R-squared:  0.7845 
## F-statistic: 18.47 on 5 and 19 DF,  p-value: 1.044e-06

model4<- class="identifier" color="#000000" font="">lm(Price~AGST+HarvestRain+WinterRain+Age,data=wine)
summary(model4)

## 
## Call:
## lm(formula = Price ~ AGST + HarvestRain + WinterRain + Age, data = wine)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.45470 -0.24273  0.00752  0.19773  0.53637 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -3.4299802  1.7658975  -1.942 0.066311 .  
## AGST         0.6072093  0.0987022   6.152  5.2e-06 ***
## HarvestRain -0.0039715  0.0008538  -4.652 0.000154 ***
## WinterRain   0.0010755  0.0005073   2.120 0.046694 *  
## Age          0.0239308  0.0080969   2.956 0.007819 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.295 on 20 degrees of freedom
## Multiple R-squared:  0.8286, Adjusted R-squared:  0.7943 
## F-statistic: 24.17 on 4 and 20 DF,  p-value: 2.036e-07

cor(wine$WinterRain,wine$Price)

## [1] 0.1366505

cor(wine$Age,wine$FrancePop)

## [1] -0.9944851

cor(wine)

##                    Year      Price   WinterRain        AGST HarvestRain
## Year         1.00000000 -0.4477679  0.016970024 -0.24691585  0.02800907
## Price       -0.44776786  1.0000000  0.136650547  0.65956286 -0.56332190
## WinterRain   0.01697002  0.1366505  1.000000000 -0.32109061 -0.27544085
## AGST        -0.24691585  0.6595629 -0.321090611  1.00000000 -0.06449593
## HarvestRain  0.02800907 -0.5633219 -0.275440854 -0.06449593  1.00000000
## Age         -1.00000000  0.4477679 -0.016970024  0.24691585 -0.02800907
## FrancePop    0.99448510 -0.4668616 -0.001621627 -0.25916227  0.04126439
##                     Age    FrancePop
## Year        -1.00000000  0.994485097
## Price        0.44776786 -0.466861641
## WinterRain  -0.01697002 -0.001621627
## AGST         0.24691585 -0.259162274
## HarvestRain -0.02800907  0.041264394
## Age          1.00000000 -0.994485097
## FrancePop   -0.99448510  1.000000000

model5<-lm>(Price~AGST+HarvestRain+WinterRain,data=wine)
summary(model5)

## 
## Call:
## lm(formula = Price ~ AGST + HarvestRain + WinterRain, data = wine)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -0.67472 -0.12958  0.01973  0.20751  0.63846 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) -4.3016263  2.0366743  -2.112 0.046831 *  
## AGST         0.6810242  0.1117011   6.097 4.75e-06 ***
## HarvestRain -0.0039481  0.0009987  -3.953 0.000726 ***
## WinterRain   0.0011765  0.0005920   1.987 0.060097 .  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 0.345 on 21 degrees of freedom
## Multiple R-squared:  0.7537, Adjusted R-squared:  0.7185 
## F-statistic: 21.42 on 3 and 21 DF,  p-value: 1.359e-06

wineTest<- class="identifier" color="#000000" font="">read.csv('Wine_test.csv')
str(wineTest)

## 'data.frame':    2 obs. of  7 variables:
##  $ Year       : int  1979 1980
##  $ Price      : num  6.95 6.5
##  $ WinterRain : int  717 578
##  $ AGST       : num  16.2 16
##  $ HarvestRain: int  122 74
##  $ Age        : int  4 3
##  $ FrancePop  : num  54836 55110

predictTest<- span=""> predict(model4,newdata=wineTest)
predictTest

##        1        2 
## 6.768925 6.684910

SST<- span=""> sum((wineTest$Price-mean(wine$Price))^2)
1-SSE/SST

## [1] -4.140916

Model Selection

2014-11-25T04:31:00.001+09:00

library(ISLR)

## Warning: package 'ISLR' was built under R version 3.1.2

summary(Hitters)

##      AtBat            Hits         HmRun            Runs       
##  Min.   : 16.0   Min.   :  1   Min.   : 0.00   Min.   :  0.00  
##  1st Qu.:255.2   1st Qu.: 64   1st Qu.: 4.00   1st Qu.: 30.25  
##  Median :379.5   Median : 96   Median : 8.00   Median : 48.00  
##  Mean   :380.9   Mean   :101   Mean   :10.77   Mean   : 50.91  
##  3rd Qu.:512.0   3rd Qu.:137   3rd Qu.:16.00   3rd Qu.: 69.00  
##  Max.   :687.0   Max.   :238   Max.   :40.00   Max.   :130.00  
##                                                                
##       RBI             Walks            Years            CAtBat       
##  Min.   :  0.00   Min.   :  0.00   Min.   : 1.000   Min.   :   19.0  
##  1st Qu.: 28.00   1st Qu.: 22.00   1st Qu.: 4.000   1st Qu.:  816.8  
##  Median : 44.00   Median : 35.00   Median : 6.000   Median : 1928.0  
##  Mean   : 48.03   Mean   : 38.74   Mean   : 7.444   Mean   : 2648.7  
##  3rd Qu.: 64.75   3rd Qu.: 53.00   3rd Qu.:11.000   3rd Qu.: 3924.2  
##  Max.   :121.00   Max.   :105.00   Max.   :24.000   Max.   :14053.0  
##                                                                      
##      CHits            CHmRun           CRuns             CRBI        
##  Min.   :   4.0   Min.   :  0.00   Min.   :   1.0   Min.   :   0.00  
##  1st Qu.: 209.0   1st Qu.: 14.00   1st Qu.: 100.2   1st Qu.:  88.75  
##  Median : 508.0   Median : 37.50   Median : 247.0   Median : 220.50  
##  Mean   : 717.6   Mean   : 69.49   Mean   : 358.8   Mean   : 330.12  
##  3rd Qu.:1059.2   3rd Qu.: 90.00   3rd Qu.: 526.2   3rd Qu.: 426.25  
##  Max.   :4256.0   Max.   :548.00   Max.   :2165.0   Max.   :1659.00  
##                                                                      
##      CWalks        League  Division    PutOuts          Assists     
##  Min.   :   0.00   A:175   E:157    Min.   :   0.0   Min.   :  0.0  
##  1st Qu.:  67.25   N:147   W:165    1st Qu.: 109.2   1st Qu.:  7.0  
##  Median : 170.50                    Median : 212.0   Median : 39.5  
##  Mean   : 260.24                    Mean   : 288.9   Mean   :106.9  
##  3rd Qu.: 339.25                    3rd Qu.: 325.0   3rd Qu.:166.0  
##  Max.   :1566.00                    Max.   :1378.0   Max.   :492.0  
##                                                                     
##      Errors          Salary       NewLeague
##  Min.   : 0.00   Min.   :  67.5   A:176    
##  1st Qu.: 3.00   1st Qu.: 190.0   N:146    
##  Median : 6.00   Median : 425.0            
##  Mean   : 8.04   Mean   : 535.9            
##  3rd Qu.:11.00   3rd Qu.: 750.0            
##  Max.   :32.00   Max.   :2460.0            
##                  NA's   :59

There are some missing values here,so before we proceed we will remove them:

Hitters<- span=""> na.omit(Hitters)
with(Hitters,sum(is.na(Salary)))

## [1] 0

Best Subset regression

We will now use the package leaps to evaluate all the best-subset models.

library(leaps)
regfit.full<- span=""> regsubsets(Salary~.,data=Hitters)
summary(regfit.full)

## Subset selection object
## Call: regsubsets.formula(Salary ~ ., data = Hitters)
## 19 Variables  (and intercept)
##            Forced in Forced out
## AtBat          FALSE      FALSE
## Hits           FALSE      FALSE
## HmRun          FALSE      FALSE
## Runs           FALSE      FALSE
## RBI            FALSE      FALSE
## Walks          FALSE      FALSE
## Years          FALSE      FALSE
## CAtBat         FALSE      FALSE
## CHits          FALSE      FALSE
## CHmRun         FALSE      FALSE
## CRuns          FALSE      FALSE
## CRBI           FALSE      FALSE
## CWalks         FALSE      FALSE
## LeagueN        FALSE      FALSE
## DivisionW      FALSE      FALSE
## PutOuts        FALSE      FALSE
## Assists        FALSE      FALSE
## Errors         FALSE      FALSE
## NewLeagueN     FALSE      FALSE
## 1 subsets of each size up to 8
## Selection Algorithm: exhaustive
##          AtBat Hits HmRun Runs RBI Walks Years CAtBat CHits CHmRun CRuns
## 1  ( 1 ) " "   " "  " "   " "  " " " "   " "   " "    " "   " "    " "  
## 2  ( 1 ) " "   "*"  " "   " "  " " " "   " "   " "    " "   " "    " "  
## 3  ( 1 ) " "   "*"  " "   " "  " " " "   " "   " "    " "   " "    " "  
## 4  ( 1 ) " "   "*"  " "   " "  " " " "   " "   " "    " "   " "    " "  
## 5  ( 1 ) "*"   "*"  " "   " "  " " " "   " "   " "    " "   " "    " "  
## 6  ( 1 ) "*"   "*"  " "   " "  " " "*"   " "   " "    " "   " "    " "  
## 7  ( 1 ) " "   "*"  " "   " "  " " "*"   " "   "*"    "*"   "*"    " "  
## 8  ( 1 ) "*"   "*"  " "   " "  " " "*"   " "   " "    " "   "*"    "*"  
##          CRBI CWalks LeagueN DivisionW PutOuts Assists Errors NewLeagueN
## 1  ( 1 ) "*"  " "    " "     " "       " "     " "     " "    " "       
## 2  ( 1 ) "*"  " "    " "     " "       " "     " "     " "    " "       
## 3  ( 1 ) "*"  " "    " "     " "       "*"     " "     " "    " "       
## 4  ( 1 ) "*"  " "    " "     "*"       "*"     " "     " "    " "       
## 5  ( 1 ) "*"  " "    " "     "*"       "*"     " "     " "    " "       
## 6  ( 1 ) "*"  " "    " "     "*"       "*"     " "     " "    " "       
## 7  ( 1 ) " "  " "    " "     "*"       "*"     " "     " "    " "       
## 8  ( 1 ) " "  "*"    " "     "*"       "*"     " "     " "    " "

It gives by default best-subsets up to size 8; lets increase that to 19, i.e. all the variables

regfit.full<- span=""> regsubsets(Salary~.,data=Hitters,nvmax=19)
reg.summary<- span=""> summary(regfit.full)
names(reg.summary)

## [1] "which"  "rsq"    "rss"    "adjr2"  "cp"     "bic"    "outmat" "obj"

plot(reg.summary$cp,xlab='Number of Variables',ylab='Cp')
which.min(reg.summary$cp)

## [1] 10

points(10,reg.summary$cp[10],pch=20,col='red')

There is a plot method for the regsubsetsobject

plot(regfit.full,scale='Cp')

coef(regfit.full,10)

##  (Intercept)        AtBat         Hits        Walks       CAtBat 
##  162.5354420   -2.1686501    6.9180175    5.7732246   -0.1300798 
##        CRuns         CRBI       CWalks    DivisionW      PutOuts 
##    1.4082490    0.7743122   -0.8308264 -112.3800575    0.2973726 
##      Assists 
##    0.2831680

Resampling Method

2014-11-25T03:35:00.001+09:00

require(ISLR)

## Loading required package: ISLR

## Warning: package 'ISLR' was built under R version 3.1.2

#install.packages('boot')
require(boot)

## Loading required package: boot

## Warning: package 'boot' was built under R version 3.1.2

?cv.glm

## starting httpd help server ... done

plot(mpg~horsepower,data=Auto)

##LOOGV
glm.fit<- span=""> glm(mpg~horsepower,data=Auto)
cv.glm(Auto,glm.fit)$delta #pretty slow (doesnt use formula (5.2) on page 180)

## [1] 24.23151 24.23114

##Lets write a simple function to use formula (5.2)
loocv<- span=""> function(fit){
h<- span=""> lm.influence(fit)$h
mean((residuals(fit)/(1-h))^2)
}
## Now we try it out
loocv(glm.fit)

## [1] 24.23151

cv.error<- span=""> rep(0,5)
degree<- span=""> 1:5
for(d in degree){
glm.fit<- span=""> glm(mpg~poly(horsepower,d),data=Auto)
cv.error[d]<- span=""> loocv(glm.fit)
}
plot(degree,cv.error,type='b')
## 10-fold CV
cv.error10<- span=""> rep(0,5)
for(d in degree){
glm.fit<- span=""> glm(mpg~poly(horsepower,d),data=Auto)
cv.error10[d]<- span=""> cv.glm(Auto,glm.fit,K=10)$delta[1]
}
lines(degree,cv.error10,type='b',col='red')

## Bootstrap
## Minimum risk investment - Section 5.2
alpha<- span=""> function(x,y){
vx<- span=""> var(x)
vy<- span=""> var(y)
cxy<- span=""> cov(x,y)
(vy-cxy)/(vx+vy-2+cxy)
}
alpha(Portfolio$X,Portfolio$Y)

## [1] 0.6413231

## What is the standard error of alpha?
alpha.fn<- span=""> function(data,index){
with(data[index,],alpha(X,Y))
}
alpha.fn(Portfolio,1:100)

## [1] 0.6413231

set.seed(1)
alpha.fn(Portfolio,sample(1:100,100,replace=T))

## [1] 0.6727862

boot.out<- span=""> boot(Portfolio,alpha.fn,R=1000)
boot.out

## 
## ORDINARY NONPARAMETRIC BOOTSTRAP
## 
## 
## Call:
## boot(data = Portfolio, statistic = alpha.fn, R = 1000)
## 
## 
## Bootstrap Statistics :
##      original     bias    std. error
## t1* 0.6413231 0.08995096   0.4226578

plot(boot.out)

Classfication

2014-11-25T02:53:00.000+09:00

require(ISLR)

## Loading required package: ISLR

## Warning: package 'ISLR' was built under R version 3.1.2

names(Smarket)

## [1] "Year"      "Lag1"      "Lag2"      "Lag3"      "Lag4"      "Lag5"     
## [7] "Volume"    "Today"     "Direction"

summary(Smarket)

##       Year           Lag1                Lag2          
##  Min.   :2001   Min.   :-4.922000   Min.   :-4.922000  
##  1st Qu.:2002   1st Qu.:-0.639500   1st Qu.:-0.639500  
##  Median :2003   Median : 0.039000   Median : 0.039000  
##  Mean   :2003   Mean   : 0.003834   Mean   : 0.003919  
##  3rd Qu.:2004   3rd Qu.: 0.596750   3rd Qu.: 0.596750  
##  Max.   :2005   Max.   : 5.733000   Max.   : 5.733000  
##       Lag3                Lag4                Lag5         
##  Min.   :-4.922000   Min.   :-4.922000   Min.   :-4.92200  
##  1st Qu.:-0.640000   1st Qu.:-0.640000   1st Qu.:-0.64000  
##  Median : 0.038500   Median : 0.038500   Median : 0.03850  
##  Mean   : 0.001716   Mean   : 0.001636   Mean   : 0.00561  
##  3rd Qu.: 0.596750   3rd Qu.: 0.596750   3rd Qu.: 0.59700  
##  Max.   : 5.733000   Max.   : 5.733000   Max.   : 5.73300  
##      Volume           Today           Direction 
##  Min.   :0.3561   Min.   :-4.922000   Down:602  
##  1st Qu.:1.2574   1st Qu.:-0.639500   Up  :648  
##  Median :1.4229   Median : 0.038500             
##  Mean   :1.4783   Mean   : 0.003138             
##  3rd Qu.:1.6417   3rd Qu.: 0.596750             
##  Max.   :3.1525   Max.   : 5.733000

?Smarket

## starting httpd help server ... done

pairs(Smarket,col=Smarket$Direction)

#Logical regression
glm.fit<- span=""> glm(Direction~Lag1+Lag2+Lag3+Lag4+Lag5+Volume,data=Smarket,family=binomial)
summary(glm.fit)

## 
## Call:
## glm(formula = Direction ~ Lag1 + Lag2 + Lag3 + Lag4 + Lag5 + 
##     Volume, family = binomial, data = Smarket)
## 
## Deviance Residuals: 
##    Min      1Q  Median      3Q     Max  
## -1.446  -1.203   1.065   1.145   1.326  
## 
## Coefficients:
##              Estimate Std. Error z value Pr(>|z|)
## (Intercept) -0.126000   0.240736  -0.523    0.601
## Lag1        -0.073074   0.050167  -1.457    0.145
## Lag2        -0.042301   0.050086  -0.845    0.398
## Lag3         0.011085   0.049939   0.222    0.824
## Lag4         0.009359   0.049974   0.187    0.851
## Lag5         0.010313   0.049511   0.208    0.835
## Volume       0.135441   0.158360   0.855    0.392
## 
## (Dispersion parameter for binomial family taken to be 1)
## 
##     Null deviance: 1731.2  on 1249  degrees of freedom
## Residual deviance: 1727.6  on 1243  degrees of freedom
## AIC: 1741.6
## 
## Number of Fisher Scoring iterations: 3

glm.probs<- span=""> predict(glm.fit,type='response')
summary(glm.probs)

##    Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
##  0.4084  0.5020  0.5180  0.5184  0.5338  0.6486

glm.probs[1:5]

##         1         2         3         4         5 
## 0.5070841 0.4814679 0.4811388 0.5152224 0.5107812

attach(Smarket)
glm.pred<- span=""> ifelse(glm.probs>.5,'Up','Down')
table(glm.pred,Direction)

##         Direction
## glm.pred Down  Up
##     Down  145 141
##     Up    457 507

mean(glm.pred==Direction)

## [1] 0.5216

# Make training and test set
train<- span=""> Year<2005
glm.fit=glm(Direction~Lag1+Lag2+Lag3+Lag4+Lag5+Volume,
data=Smarket,family=binomial,subset=train)
glm.probs<- span=""> predict(glm.fit,newdata=Smarket[!train,],type='response')
glm.pred<- span=""> ifelse(glm.probs>.5,'Up','Down')
Direction.2005<- span=""> Smarket$Direction[!train]
table(glm.pred,Direction.2005)

##         Direction.2005
## glm.pred Down Up
##     Down   77 97
##     Up     34 44

mean(glm.pred==Direction.2005)

## [1] 0.4801587

#Fit smaller model
glm.git<- span=""> glm(Direction~Lag1+Lag2,data=Smarket,family=binomial,subset=train)
glm.probs<- span=""> predict(glm.fit,newdata=Smarket[!train,],type='response')
glm.pred<- span=""> ifelse(glm.probs>.5,"Up","Down")
table(glm.pred,Direction.2005)

##         Direction.2005
## glm.pred Down Up
##     Down   77 97
##     Up     34 44

mean(glm.pred==Direction.2005)

## [1] 0.4801587

106/(76+106)

## [1] 0.5824176

Linear Reggression

2014-11-24T04:33:00.003+09:00

# Chapter 3 Lab: Linear Regression
library(MASS)
library(ISLR)

## Warning: package 'ISLR' was built under R version 3.1.2

names(Boston)

##  [1] "crim"    "zn"      "indus"   "chas"    "nox"     "rm"      "age"    
##  [8] "dis"     "rad"     "tax"     "ptratio" "black"   "lstat"   "medv"

?Boston

## starting httpd help server ... done

plot(medv~lstat,data=Boston)
fit1<- span=""> lm(medv~lstat,data=Boston)
fit1

## 
## Call:
## lm(formula = medv ~ lstat, data = Boston)
## 
## Coefficients:
## (Intercept)        lstat  
##       34.55        -0.95

summary(fit1)

## 
## Call:
## lm(formula = medv ~ lstat, data = Boston)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -15.168  -3.990  -1.318   2.034  24.500 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 34.55384    0.56263   61.41   <2e-16 ---="" -0.95005="" -24.53="" 0.001="" 0.01="" 0.03873="" 0.05="" 0.1="" 0.5432="" 0.5441="" 0="" 1="" 2.2e-16="" 504="" 6.216="" 601.6="" adjusted="" and="" code="" codes:="" degrees="" df="" e-16="" error:="" f-statistic:="" freedom="" lstat="" multiple="" of="" on="" p-value:="" r-squared:="" residual="" signif.="" standard="">

abline(fit1,col='red')

names(fit1)

##  [1] "coefficients"  "residuals"     "effects"       "rank"         
##  [5] "fitted.values" "assign"        "qr"            "df.residual"  
##  [9] "xlevels"       "call"          "terms"         "model"

confint(fit1)

##                 2.5 %     97.5 %
## (Intercept) 33.448457 35.6592247
## lstat       -1.026148 -0.8739505

predict(fit1,data.frame(lstat=c(5,10,15)),interval='confidence')

##        fit      lwr      upr
## 1 29.80359 29.00741 30.59978
## 2 25.05335 24.47413 25.63256
## 3 20.30310 19.73159 20.87461

##Multiple linear regression
fit2<- span=""> lm(medv~lstat+age,data=Boston)
summary(fit2)

## 
## Call:
## lm(formula = medv ~ lstat + age, data = Boston)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -15.981  -3.978  -1.283   1.968  23.158 
## 
## Coefficients:
##             Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 33.22276    0.73085  45.458  < 2e-16 ***
## lstat       -1.03207    0.04819 -21.416  < 2e-16 ***
## age          0.03454    0.01223   2.826  0.00491 ** 
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.173 on 503 degrees of freedom
## Multiple R-squared:  0.5513, Adjusted R-squared:  0.5495 
## F-statistic:   309 on 2 and 503 DF,  p-value: < 2.2e-16

fit3<- span=""> lm(medv~.,Boston)
summary(fit3)

## 
## Call:
## lm(formula = medv ~ ., data = Boston)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -15.595  -2.730  -0.518   1.777  26.199 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  3.646e+01  5.103e+00   7.144 3.28e-12 ***
## crim        -1.080e-01  3.286e-02  -3.287 0.001087 ** 
## zn           4.642e-02  1.373e-02   3.382 0.000778 ***
## indus        2.056e-02  6.150e-02   0.334 0.738288    
## chas         2.687e+00  8.616e-01   3.118 0.001925 ** 
## nox         -1.777e+01  3.820e+00  -4.651 4.25e-06 ***
## rm           3.810e+00  4.179e-01   9.116  < 2e-16 ***
## age          6.922e-04  1.321e-02   0.052 0.958229    
## dis         -1.476e+00  1.995e-01  -7.398 6.01e-13 ***
## rad          3.060e-01  6.635e-02   4.613 5.07e-06 ***
## tax         -1.233e-02  3.760e-03  -3.280 0.001112 ** 
## ptratio     -9.527e-01  1.308e-01  -7.283 1.31e-12 ***
## black        9.312e-03  2.686e-03   3.467 0.000573 ***
## lstat       -5.248e-01  5.072e-02 -10.347  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.745 on 492 degrees of freedom
## Multiple R-squared:  0.7406, Adjusted R-squared:  0.7338 
## F-statistic: 108.1 on 13 and 492 DF,  p-value: < 2.2e-16

par(mfrow=c(2,2))
plot(fit3)

fit4<- span=""> update(fit3,~.-age-indus)
summary(fit4)

## 
## Call:
## lm(formula = medv ~ crim + zn + chas + nox + rm + dis + rad + 
##     tax + ptratio + black + lstat, data = Boston)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -15.5984  -2.7386  -0.5046   1.7273  26.2373 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept)  36.341145   5.067492   7.171 2.73e-12 ***
## crim         -0.108413   0.032779  -3.307 0.001010 ** 
## zn            0.045845   0.013523   3.390 0.000754 ***
## chas          2.718716   0.854240   3.183 0.001551 ** 
## nox         -17.376023   3.535243  -4.915 1.21e-06 ***
## rm            3.801579   0.406316   9.356  < 2e-16 ***
## dis          -1.492711   0.185731  -8.037 6.84e-15 ***
## rad           0.299608   0.063402   4.726 3.00e-06 ***
## tax          -0.011778   0.003372  -3.493 0.000521 ***
## ptratio      -0.946525   0.129066  -7.334 9.24e-13 ***
## black         0.009291   0.002674   3.475 0.000557 ***
## lstat        -0.522553   0.047424 -11.019  < 2e-16 ***
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 4.736 on 494 degrees of freedom
## Multiple R-squared:  0.7406, Adjusted R-squared:  0.7348 
## F-statistic: 128.2 on 11 and 494 DF,  p-value: < 2.2e-16

###Nonlinear terms and Interactons
fit5<- span=""> lm(medv~lstat*age,Boston)
summary(fit5)

## 
## Call:
## lm(formula = medv ~ lstat * age, data = Boston)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -15.806  -4.045  -1.333   2.085  27.552 
## 
## Coefficients:
##               Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 36.0885359  1.4698355  24.553  < 2e-16 ***
## lstat       -1.3921168  0.1674555  -8.313 8.78e-16 ***
## age         -0.0007209  0.0198792  -0.036   0.9711    
## lstat:age    0.0041560  0.0018518   2.244   0.0252 *  
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 6.149 on 502 degrees of freedom
## Multiple R-squared:  0.5557, Adjusted R-squared:  0.5531 
## F-statistic: 209.3 on 3 and 502 DF,  p-value: < 2.2e-16

fit6<- span=""> lm(medv~lstat+I(lstat^2),Boston)
summary(fit6)

## 
## Call:
## lm(formula = medv ~ lstat + I(lstat^2), data = Boston)
## 
## Residuals:
##      Min       1Q   Median       3Q      Max 
## -15.2834  -3.8313  -0.5295   2.3095  25.4148 
## 
## Coefficients:
##              Estimate Std. Error t value Pr(>|t|)    
## (Intercept) 42.862007   0.872084   49.15   <2e-16 ---="" -18.84="" -2.332821="" 0.001="" 0.003745="" 0.01="" 0.043547="" 0.05="" 0.123803="" 0.1="" 0.6393="" 0.6407="" 0="" 11.63="" 1="" 2.2e-16="" 2="" 448.5="" 5.524="" 503="" adjusted="" and="" code="" codes:="" degrees="" df="" e-16="" error:="" f-statistic:="" freedom="" i="" lstat="" multiple="" of="" on="" p-value:="" r-squared:="" residual="" signif.="" standard="">

attach(Boston)
par(mfrow=c(1,1))
plot(medv~lstat)
points(lstat,fitted(fit6),col='red',pch=20)
fit7<- span=""> lm(medv~poly(lstat,4))
points(lstat,fitted(fit7),col='blue',pch=20)

plot(1:20,1:20,pch=1:20,cex=2)

###Qualitative predictors
fix(Carseats)
names(Carseats)

##  [1] "Sales"       "CompPrice"   "Income"      "Advertising" "Population" 
##  [6] "Price"       "ShelveLoc"   "Age"         "Education"   "Urban"      
## [11] "US"

summary(Carseats)

##      Sales          CompPrice       Income        Advertising    
##  Min.   : 0.000   Min.   : 77   Min.   : 21.00   Min.   : 0.000  
##  1st Qu.: 5.390   1st Qu.:115   1st Qu.: 42.75   1st Qu.: 0.000  
##  Median : 7.490   Median :125   Median : 69.00   Median : 5.000  
##  Mean   : 7.496   Mean   :125   Mean   : 68.66   Mean   : 6.635  
##  3rd Qu.: 9.320   3rd Qu.:135   3rd Qu.: 91.00   3rd Qu.:12.000  
##  Max.   :16.270   Max.   :175   Max.   :120.00   Max.   :29.000  
##    Population        Price        ShelveLoc        Age       
##  Min.   : 10.0   Min.   : 24.0   Bad   : 96   Min.   :25.00  
##  1st Qu.:139.0   1st Qu.:100.0   Good  : 85   1st Qu.:39.75  
##  Median :272.0   Median :117.0   Medium:219   Median :54.50  
##  Mean   :264.8   Mean   :115.8                Mean   :53.32  
##  3rd Qu.:398.5   3rd Qu.:131.0                3rd Qu.:66.00  
##  Max.   :509.0   Max.   :191.0                Max.   :80.00  
##    Education    Urban       US     
##  Min.   :10.0   No :118   No :142  
##  1st Qu.:12.0   Yes:282   Yes:258  
##  Median :14.0                      
##  Mean   :13.9                      
##  3rd Qu.:16.0                      
##  Max.   :18.0

fit1<- span=""> lm(Sales~.+Income:Advertising+Age:Price,Carseats)
summary(fit1)

## 
## Call:
## lm(formula = Sales ~ . + Income:Advertising + Age:Price, data = Carseats)
## 
## Residuals:
##     Min      1Q  Median      3Q     Max 
## -2.9208 -0.7503  0.0177  0.6754  3.3413 
## 
## Coefficients:
##                      Estimate Std. Error t value Pr(>|t|)    
## (Intercept)         6.5755654  1.0087470   6.519 2.22e-10 ***
## CompPrice           0.0929371  0.0041183  22.567  < 2e-16 ***
## Income              0.0108940  0.0026044   4.183 3.57e-05 ***
## Advertising         0.0702462  0.0226091   3.107 0.002030 ** 
## Population          0.0001592  0.0003679   0.433 0.665330    
## Price              -0.1008064  0.0074399 -13.549  < 2e-16 ***
## ShelveLocGood       4.8486762  0.1528378  31.724  < 2e-16 ***
## ShelveLocMedium     1.9532620  0.1257682  15.531  < 2e-16 ***
## Age                -0.0579466  0.0159506  -3.633 0.000318 ***
## Education          -0.0208525  0.0196131  -1.063 0.288361    
## UrbanYes            0.1401597  0.1124019   1.247 0.213171    
## USYes              -0.1575571  0.1489234  -1.058 0.290729    
## Income:Advertising  0.0007510  0.0002784   2.698 0.007290 ** 
## Price:Age           0.0001068  0.0001333   0.801 0.423812    
## ---
## Signif. codes:  0 '***' 0.001 '**' 0.01 '*' 0.05 '.' 0.1 ' ' 1
## 
## Residual standard error: 1.011 on 386 degrees of freedom
## Multiple R-squared:  0.8761, Adjusted R-squared:  0.8719 
## F-statistic:   210 on 13 and 386 DF,  p-value: < 2.2e-16

contrasts(Carseats$ShelveLoc)

##        Good Medium
## Bad       0      0
## Good      1      0
## Medium    0      1

###Writing R function
regplot<- span=""> function(x,y){
fit<- span=""> lm(y~x)
plot(x,y)
abline(fit,col='red')
}
attach(Carseats)
regplot(Price,Sales)

regplot<- span=""> function(x,y,...){
fit<- span=""> lm(y~x)
plot(x,y,...)
abline(fit,col='red')
}
regplot(Price,Sales,xlab='Price',ylab='Sales',col='blue',pch=20)

Alfa Investor In A Good Mood

Nikkei225

Nikkei 225

nikkei 225

Return Calculus R

4銘柄のボラティリティ

`sd(NipponSheetGlassCompany_TSE5202)`

`## [1] 20.12854`

`sd(PolaOrbisHoldings_TSE4927)`

`## [1] 209.8629`

`sd(YamahaMotor_TSE7272)`

`## [1] 187.3384`

`sd(TakedaPharmaceuticalCompany_TSE4502)`

`## [1] 571.4873`

μ、σ、histgram

`par(mfcol=c(2,2)) hist(NipponSheetGlassCompany_TSE5202,main="NipponSheetGlassCompany_TSE5202") hist(PolaOrbisHoldings_TSE4927,main="PolaOrbisHoldings_TSE4927") hist(YamahaMotor_TSE7272,main="YamahaMotor_TSE7272") hist(TakedaPharmaceuticalCompany_TSE4502,main="TakedaPharmaceuticalCompany_TSE4502")`

散布図行列

Shapiro-Wilk検定

連の検定(runs test)

analyze four historical time series

The function normality_tests combines three different statistical tests:

the German DAX index :^GDAXI the American S&P 500 index:

^GSPC YHOO: Yahoo MSFT: Micro Soft

Index start at 100

Installed pymc library Python

pymc 2.3.4 --Python Library(Markov Chain Monte Carlo sampling toolkit)

Go to the download site here .

Time dependency

Time dependency

Wine Prediction

Wine_Regression

Model Selection

Model Selection

Best Subset regression

Resampling Method

Resampling Method

Classfication

Classfication.Rmd

Linear Reggression

	^GDAXI	^GSPC	YHOO	MSFT
Date
2006-01-03	NaN	NaN	NaN	NaN
2006-01-04	0.011460	0.003666	0.001466	0.004832
2006-01-05	-0.001284	0.000016	0.013576	0.000741
2006-01-06	0.003581	0.009356	0.039656	-0.002968
2006-01-09	0.000143	0.003650	0.004848	-0.001860