2014年9月12日から過去三カ月間の日経平均株価終値とドル円公示レートの関係(データ出所:株式会社日本経済新聞社、株式会社みずほ銀行)

nikkei_stock_average_daily_en.csv_NIKKEIClose quote.csv_USD

CCF_NIKKEIClose-USD


File:nikkei_stock_average_daily_en.csv ColumnName:NIKKEIClose Period:2014/6/13-2014/9/12
1-Unit Root Test H0:Contains a Unit Root(Random Walk) 
- Original p-value=0.3432429 Lag=0
- 1st Difference p-value=0.01 Lag=0

2-Summary 
- Original 
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  14780   15280   15370   15390   15530   15950 
- 1st Difference 
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
-454.00  -47.41   16.73   13.29   72.95  352.10 

3-Other Descriptive Statistics 
- Original 
-- n= 65  Unbiased Variance= 49941.3  Standard Variance= 223.4755
- 1st Difference 
-- n= 64  Unbiased Variance= 13740.03  Standard Variance= 117.2179
 
4-ARIMA 
- Original
 ARIMA(2,1,2) with drift         : 1e+20
 ARIMA(0,1,0) with drift         : 782.0625
 ARIMA(1,1,0) with drift         : 782.7037
 ARIMA(0,1,1) with drift         : 782.9507
 ARIMA(1,1,1) with drift         : 784.3565
 ARIMA(0,1,0)                    : 793.2428

 Best model: ARIMA(0,1,0) with drift         

Series: dataset[, ccc] 
ARIMA(0,1,0) with drift         

Coefficients:
        drift
      13.2883
s.e.  14.6522

sigma^2 estimated as 13740:  log likelihood=-389.03
AIC=782.06   AICc=782.26   BIC=786.38

5-Normality Test for Residual of ARIMA H0:Obey a Normal Distribution 
- Original
	Shapiro-Wilk normality test

data:  Result.arima$res
W = 0.9448, p-value = 0.005868

6-Part of Original 
          date NIKKEIClose
847 2014-06-13    15097.84
848 2014-06-16    14933.29
849 2014-06-17    14975.97
850 2014-06-18    15115.80
851 2014-06-19    15361.16
852 2014-06-20    15349.42
853 2014-06-23    15369.28
854 2014-06-24    15376.24
855 2014-06-25    15266.61
856 2014-06-26    15308.49

 





File:quote.csv ColumnName:USD Period:2014/6/13-2014/9/12
1-Unit Root Test H0:Contains a Unit Root(Random Walk) 
- Original p-value=0.99 Lag=0
- 1st Difference p-value=0.01 Lag=0

2-Summary 
- Original 
   Min. 1st Qu.  Median    Mean 3rd Qu.    Max. 
  101.3   101.7   102.2   102.7   103.7   107.2 
- 1st Difference 
    Min.  1st Qu.   Median     Mean  3rd Qu.     Max. 
-0.46000 -0.09500  0.04500  0.08469  0.24000  0.93000 

3-Other Descriptive Statistics 
- Original 
-- n= 65  Unbiased Variance= 2.121764  Standard Variance= 1.456628
- 1st Difference 
-- n= 64  Unbiased Variance= 0.09465704  Standard Variance= 0.3076638
 
4-ARIMA 
- Original
 ARIMA(2,2,2)                    : 34.0331
 ARIMA(0,2,0)                    : 76.18334
 ARIMA(1,2,0)                    : 41.13219
 ARIMA(0,2,1)                    : 35.40463
 ARIMA(1,2,2)                    : 35.22105
 ARIMA(3,2,2)                    : 32.20557
 ARIMA(3,2,1)                    : 32.60964
 ARIMA(3,2,3)                    : 1e+20
 ARIMA(2,2,1)                    : 34.18205
 ARIMA(4,2,3)                    : 35.39768
 ARIMA(3,2,2)                    : 32.20557
 ARIMA(4,2,2)                    : 36.50742

 Best model: ARIMA(3,2,2)                    

Series: dataset[, ccc] 
ARIMA(3,2,2)                    

Coefficients:
         ar1     ar2      ar3      ma1     ma2
      0.8113  0.2700  -0.4034  -1.8549  0.8954
s.e.  0.1412  0.1475   0.1227   0.1133  0.1122

sigma^2 estimated as 0.07665:  log likelihood=-10.1
AIC=32.21   AICc=33.71   BIC=45.06

5-Normality Test for Residual of ARIMA H0:Obey a Normal Distribution 
- Original
	Shapiro-Wilk normality test

data:  Result.arima$res
W = 0.9804, p-value = 0.392

6-Part of Original 
           date    USD
2997 2014-06-13 101.80
2998 2014-06-16 101.98
2999 2014-06-17 101.95
3000 2014-06-18 102.18
3001 2014-06-19 102.01
3002 2014-06-20 101.90
3003 2014-06-23 102.10
3004 2014-06-24 101.83
3005 2014-06-25 101.96
3006 2014-06-26 101.70

 



co-integration test between NIKKEIClose and USD H0:no cointegrated.error term is not stationary
######################################## 
# Phillips and Ouliaris Unit Root Test # 
######################################## 

Test of type Pu 
detrending of series with constant and linear trend 


Call:
lm(formula = z[, 1] ~ z[, -1] + trd)

Residuals:
     Min       1Q   Median       3Q      Max 
-1.43004 -0.53166 -0.07938  0.38567  1.75644 

Coefficients:
             Estimate Std. Error t value Pr(>|t|)    
(Intercept) 6.092e+01  8.242e+00   7.391 4.51e-10 ***
z[, -1]     2.626e-03  5.443e-04   4.826 9.47e-06 ***
trd         4.227e-02  6.433e-03   6.570 1.19e-08 ***
---
Signif. codes:  0 ‘***’ 0.001 ‘**’ 0.01 ‘*’ 0.05 ‘.’ 0.1 ‘ ’ 1

Residual standard error: 0.7368 on 62 degrees of freedom
Multiple R-squared:  0.7521,	Adjusted R-squared:  0.7441 
F-statistic: 94.07 on 2 and 62 DF,  p-value: < 2.2e-16


Value of test-statistic is: 5.422 

Critical values of Pu are:
                  10pct    5pct    1pct
critical values 41.2488 48.8439 65.1714


D.W test between NIKKEIClose and USD H0:the autocorrelation of the disturbances is 0
	Durbin-Watson test

data:  lm(AllData02[, c(nnn, 2)])
DW = 0.3695, p-value < 2.2e-16
alternative hypothesis: true autocorrelation is greater than 0






Warning messages:
1: In adf.test(diff(dataset[, ccc]), k = UnitRootLag) :
  p-value smaller than printed p-value
2: In xy.coords(x, y, xlabel, ylabel, log) :
  1 y value <= 0 omitted from logarithmic plot
3: In adf.test(dataset[, ccc], k = UnitRootLag) :
  p-value greater than printed p-value
4: In adf.test(diff(dataset[, ccc]), k = UnitRootLag) :
  p-value smaller than printed p-value
5: In xy.coords(x, y, xlabel, ylabel, log) :
  1 y value <= 0 omitted from logarithmic plot
> 
#initial condition
FirstDate<-as.Date("2014/6/13")
LastDate<-as.Date("2014/9/12")
period<-"days"
condition01<-0 #analysis each series 0:done,1:not done
condition02<-0 #mutual analysis 0:done,1:not done
UnitRootLag=0
SL<-5 #significant level for structure change
if(period=="days"){
DateStyle<-"%Y/%b/%d"
}else
if(period=="months"){
DateStyle<-"%Y/%b"  
}
#package
library(urca) # http://cran.r-project.org/web/packages/urca/urca.pdf
library(tseries) # http://cran.r-project.org/web/packages/tseries/tseries.pdf
library(forecast) # http://cran.r-project.org/web/packages/forecast/forecast.pdf
library(MASS) # http://cran.r-project.org/web/packages/MASS/MASS.pdf
library(lmtest) # http://cran.r-project.org/web/packages/lmtest/lmtest.pdf
library(systemfit) # http://cran.r-project.org/web/packages/systemfit/systemfit.pdf
#file path
username<-Sys.info()['user']
path01<-paste("C:/Users/",username,"/Desktop/R_Data_Read/",sep="")
path02<-paste("C:/Users/",username,"/Desktop/R_Graph/",sep="")
path03<-paste("C:/Users/",username,"/Desktop/R_Data_Write/",sep="")
#start analysis
setwd(path01)
AllData01<-data.frame(date=seq(FirstDate,LastDate,by=period),dummy=1)
AllData02<-AllData01
for(iii in 1:length(dir(path01))){
dataset<-read.table(file=dir(path01)[iii],header=TRUE,sep=",")
#date format
dataset[,1]<-as.Date(dataset[,1])
colnames(dataset)[1]<-"date"
dataset<-subset(dataset,FirstDate<=dataset[,1])
dataset<-subset(dataset,LastDate>=dataset[,1])
for(ccc in 2:ncol(dataset)){
dataset<-subset(dataset,dataset[,ccc]!=".")
dataset<-subset(dataset,is.na(dataset[,ccc])==F)
dataset[,ccc]<-as.double(as.character(dataset[,ccc]))
if(condition01==0){
#text part
#unit root test
p0<-adf.test(dataset[,ccc],k=UnitRootLag)$p.value	
p1<-adf.test(diff(dataset[,ccc]),k=UnitRootLag)$p.value	
cat("\n","\n","File:",dir(path01)[iii]," ColumnName:",colnames(dataset)[ccc],
" Period:",format(FirstDate,DateStyle),"-",format(LastDate,DateStyle),"\n",sep="")
cat("1-Unit Root Test H0:Contains a Unit Root(Random Walk)","\n")
# http://cran.r-project.org/web/packages/tseries/tseries.pdf#page=2
cat("- Original p-value=",p0," Lag=",UnitRootLag,"\n",sep="")
cat("- 1st Difference p-value=",p1," Lag=",UnitRootLag,"\n",sep="")
cat("\n")
cat("2-Summary","\n")
cat("- Original","\n")
print(summary(dataset[,ccc]))
cat("- 1st Difference","\n")
print(summary(diff(dataset[,ccc])))
cat("\n")
cat("3-Other Descriptive Statistics","\n")
cat("- Original","\n")
cat("-- n=",length(dataset[,ccc])," ")
cat("Unbiased Variance=",var(dataset[,ccc])," ")
cat("Standard Variance=",sd(dataset[,ccc]))
cat("\n")
cat("- 1st Difference","\n")
cat("-- n=",length(diff(dataset[,ccc]))," ")
cat("Unbiased Variance=",var(diff(dataset[,ccc]))," ")
cat("Standard Variance=",sd(diff(dataset[,ccc])))
cat("\n","\n")
cat("4-ARIMA","\n")
cat("- Original")
Result.arima<-auto.arima(dataset[,ccc],ic="aic",trace=T,stepwise=T)
# http://cran.r-project.org/web/packages/forecast/forecast.pdf#page=11
print(Result.arima)
cat("\n")
cat("5-Normality Test for Residual of ARIMA H0:Obey a Normal Distribution","\n")
cat("- Original")
print(shapiro.test(Result.arima$res))
cat("6-Part of Original","\n")
print(head(dataset,10))
cat("\n","\n")

#graph part
ymax<-max(dataset[,ccc])
ymin<-min(dataset[,ccc])
ytitle<-colnames(dataset)[ccc]
setwd(path02)
png(file=paste(dir(path01)[iii],"_",colnames(dataset)[ccc],".png",sep=""),width=1400,height=800)
par(mar=c(5,5,5,5),ps=22,cex.axis=1,cex.lab=1,cex.main=1,cex.sub=1)
layout(matrix(c(1,2,3,4,5,6,7,8,9,10,11,12),ncol=4))#or mfrow=c(3,4)
Date<-dataset[,1]
#1
plot(Date,dataset[,ccc],type="l",ylim=c(ymin,ymax),ylab=ytitle,xaxt="n",main=ytitle)
axis.Date(side=1,dataset[,1],format=DateStyle)
#2
plot(Date,dataset[,ccc]/ymax*100,type="l",ylim=c(ymin/ymax*100,100),
ylab=paste("Max=100%:",ytitle,sep=""),xaxt="n",main=paste(ytitle,"Normalization"))
axis.Date(side=1,dataset[,1],format=DateStyle)
#3
plot(dataset[,1],dataset[,ccc],xaxt="n",type="l",ylim=c(ymin,ymax),ylab=ytitle,main="Trend(caution:omit blank date)" )
abline(lm(dataset[,ccc]~dataset[,1]),col=2)
axis.Date(side=1,dataset[,1],format=DateStyle)
#4
plot(Date[-1],diff(dataset[,ccc]),xaxt="n",type="h",ylab=paste("1stDifference:",ytitle,sep=""),main="1st Difference")
axis.Date(side=1,dataset[,1][-1],format=DateStyle)
#5
#truehist(diff(dataset[,ccc]),xlab="",main=paste("1stDifference:",ytitle,sep="")) 
hist(diff(dataset[,ccc]),xlab="",main=paste("1st Difference:",ytitle,sep="")) 
#6
qqnorm(diff(dataset[,ccc]),main="1st Diff Normal Q-Q Plot")
qqline(diff(dataset[,ccc]),col=2)
#7
acf(dataset[,ccc],type="correlation",ci=c(0.9,0.95),plot=T,
    main=colnames(dataset)[ccc],lag=length(floor(dataset[,ccc]/5))) 
#8
pacf(dataset[,ccc],ci=c(0.9,0.95),plot=T,main=colnames(dataset)[ccc],lag=length(floor(dataset[,ccc]/5)))
#9
plot(forecast(Result.arima,level=c(50,95),h=floor(length(dataset[,ccc])/10)),ylab=ytitle)
#10
#Structure Change
StChg<-data.frame(id=seq(1:length(dataset[,ccc])),dummy=0,Pr=0)
for(kkk in 2:length(dataset[,ccc])){
StChg[,2]<-0
StChg[,2][kkk<=StChg[,1]]<-1
fm0<-lm(dataset[,ccc]~seq(1:length(dataset[,ccc])))
fm1<-lm(dataset[,ccc]~seq(1:length(dataset[,ccc]))+StChg[,2]+StChg[,2]*seq(1:length(dataset[,ccc])))
StChg[,3][kkk]<-anova(fm0,fm1)$Pr[2] #H0:No structure change.g = d = 0.y(i)=a+b*x(i)+g*dummy(i)+d*dummy(i)*x(i)+u(i)  
}
plot(Date,StChg[,3],type="l",main="Structure Change",xaxt="n",ylab="Pr(>F)",log="y")
axis.Date(side=1,dataset[,1],format=DateStyle)
#11
spec.pgram(dataset[,ccc],spans=c(5,5),col=2,main=paste(colnames(dataset)[ccc],"\n","Smoothed Periodogram"))
#12
spectrum(dataset[,ccc],method="ar") #method="pgram"
title(paste("\n",dir(path01)[iii]," ",colnames(dataset)[ccc]," ",
            "Period:",format(FirstDate,DateStyle),"-",format(LastDate,DateStyle),sep=""),outer=T)
dev.off()
}
}
cat("\n")
cat("\n")
cat("\n")
setwd(path01)
AllData01<-merge(AllData01,dataset,by="date",all=T,sort=T) #is.na=T
AllData02<-merge(AllData02,dataset,by="date",sort=T) #is.na=F
if(iii==length(dir(path01))){ 
AllData01<-subset(AllData01,select=-dummy)
AllData02<-subset(AllData02,select=-dummy)
if(condition02==0){
setwd(path02)
if(3<=ncol(AllData02)){
for(nnn in 3:ncol(AllData02)){
#co-integration test
cat("co-integration test between ",colnames(AllData02[2])," and ",colnames(AllData02[nnn]),
    " H0:no cointegrated.error term is not stationary",sep="")
print(summary(ca.po(AllData02[,c(nnn,2)],demean="trend")))
# http://cran.r-project.org/web/packages/urca/urca.pdf#page=14
cat("\n")
#D.W test
cat("D.W test between ",colnames(AllData02[2])," and ",colnames(AllData02[nnn]),
    " H0:the autocorrelation of the disturbances is 0",sep="")
print(dwtest(lm(AllData02[,c(nnn,2)]))) 
# http://cran.r-project.org/web/packages/lmtest/lmtest.pdf#page=13
cat("\n")
cat("\n")
cat("\n")
#CCF
png(file=paste("CCF_",colnames(AllData02[2]),"-",colnames(AllData02[nnn]),".png",sep=""),width=800,height=800)
par(mfrow=c(2,1),mar=c(5,5,5,5),ps=22,cex.axis=1,cex.lab=1,cex.main=1,cex.sub=1)
ccf(AllData02[,nnn],AllData02[,2],main=paste("CCF:",colnames(AllData02[2]),"-",colnames(AllData02[nnn])))
ccf(diff(AllData02[,nnn]),diff(AllData02[,2]),
    main=paste("CCF:1stDiff ",colnames(AllData02[2]),"-1stDiff ",colnames(AllData02[nnn])))
dev.off()
}
}
}
cat("\n")
#print(head(AllData01,4))
cat("\n")
#print(head(AllData02,4))
setwd(path03)
write.table(AllData01,file="AllData01.csv",sep=",",quote=FALSE,row.names=FALSE)
write.table(AllData02,file="AllData02.csv",sep=",",quote=FALSE,row.names=FALSE)
}
}