カオス的市場の販売予測　共立出版　大石展緒・二宮政彦先生

```r

#レンジの計算 p.46

month=c(1:12)
x=c(6,8,15,20,20,23,27,41,32,15,15,18)
Range=c()

for(i in 1:length(x)){
syouhiryou=x[1:i]
xbar=mean(syouhiryou)
sn=sd(syouhiryou)*sqrt((length(syouhiryou)-1)/length(syouhiryou))
y=cumsum(syouhiryou-xbar)
Rn=max(y)-min(y)
Range=c(Range,Rn/sn)
}

Range[1]=0

plot(month,Range,xlim=c(min(month),max(month)),ylim=c(min(Range),max(Range)),col=2,xlab="month",ylab="Range")

par(new=T)

plot(month,sqrt(month),xlim=c(min(month),max(month)),ylim=c(min(Range),max(Range)),col=3,xlab="month",ylab="Range")

#OLS

n=length(month)
t=month[3:length(month)]
Rn=Range[3:length(Range)]

HlogC=solve(matrix(c(sum(log(t)),sum(log(t)^2),n,sum(log(t))),nrow=2))%*%c(sum(log(Rn)),sum(log(t)*log(Rn)))

H=HlogC[1] #ハースト指数
C=exp(HlogC[2]) #定数
a=1/H #フラクタル次元

Rnhat=log(C)+H*log(t)
R2=cor(Rnhat,Rn)

#Rnの平均
ERn=function(N){
if(N>300){
return(sum(sqrt((N-c(1:N))/c(1:N)))/sqrt(N*pi/2))  
}else{
return(gamma((N-1)/2)*sum(sqrt((N-c(1:N))/c(1:N)))/sqrt(pi*gamma(N/2)))  
}}


```

```r

#NNモデル p.93

#表5.1

period=c(1:10)
sales_volume=c(6.3,12.6,5.7,12,6.15,11.85,5.88,12.15,5.55,12.45)
y=sales_volume
t=length(period)

D=2
dt=function(tt){
d=0
for(i in 1:D){
d=d+sum((y[t+1-i]-y[tt-i])^2)
}
return(sqrt(d))
}

dtvec=c()
for(j in (D+1):t){
dtvec=c(dtvec,dt(j))  
}

dtdat=data.frame(t=c(3:t),dtstarvec=dtvec,y=y[3:t])
dtdat=dtdat[order(dtdat$dtstarvec),]
dtstarvec=sort(dtvec)
plot(c(1:t),dtstarvec,col=2,main="予測変数確認用")
J=4 #予測変数の数
w=(1/dtstarvec[1:J])/sum(1/dtstarvec[1:J])
ythat=sum(dtdat$y[1:J]*w) #次期の予測値

```

```r

#p.119 6.1.2　相関のあるコンセンサスモデル


mu=c(61,62.5,70)
sig=sqrt(c(36,25,49))
cormat=matrix(c(1,0.6,0.5,0.6,1,0.6,0.5,0.6,1),nrow=3)
e=rep(1,length(mu))

varmat=t(cormat*sig)*sig
mustar=as.numeric(e%*%solve(varmat)%*%mu)/as.numeric(e%*%solve(varmat)%*%e)
sigstar=sqrt(1/as.numeric(e%*%solve(varmat)%*%e))

w=apply(solve(varmat),1,sum)/sum(solve(varmat))

```

```r

#p.129 RAS法の活用

Sij=matrix(c(26,20,54,20,30,41,80,12,71,34,48,23,24,69,75,14,80,38,46,25,33,104,94,24,65,26,47,24,50,41,86,15,76,33,60,26,31,77,82,15,53,32,50,25,36,50,80,15,36,30,45,23,31,54,72,12,42,33,70,23,30,54,73,12,31,24,47,26,31,74,78,11,12,20,39,15,25,37,57,10),nrow=8)

Si.=c(506,294,529,242,339,622,824,144)#製品別売上高
S.j=c(296,368,474,373,425,357,312,344,329,222)#地域別売上高
sum(Si.)==sum(S.j)

Pi.=Si./sum(Si.)
P.j=S.j/sum(S.j)
Pij=Sij/sum(Sij)

#alpha,betaを推定するための反復アルゴリズム

beta=rep(1,ncol(Sij))#初期値

ite=10
for(l in 1:ite){
beta_pre=beta
alpha=Pi./as.numeric(Pij%*%beta)  
beta=P.j/as.numeric(t(Pij)%*%alpha) 
print(sum(abs(beta_pre-beta)))
}

Pij_hat=t(t(Pij*c(alpha))*beta)
Sij_hat=Pij_hat*sum(Si.)

```

```r

#p.200 ロジスティック曲線での予想(OLS)

data("AirPassengers")
Q=c(t(AirPassengers))
t=c(1:length(Q))
t0=0


f=function(x){
value=-log(x[3]-Q)+log(x[1])-x[2]*(t-t0)
return(value)
}


Y=-log(Q)
h=0.001
lam=1
ite=10^5
sita=c(5,5,max(Q)+1)
eta=8*10^(-1)

for(l in 1:ite){
sita_pre=sita
Z=array(0,dim=c(length(Y),length(sita)))
vec=sita
for(j in 1:ncol(Z)){
vec_sub=vec;vec_sub[j]=vec_sub[j]+h
Z_sub=(f(vec_sub)-f(vec))/h
Z[,j]=Z_sub
}
beta=solve(t(Z)%*%Z+diag(lam,ncol(Z)))%*%t(Z)%*%(Y-f(sita))
sita=c(sita+eta*beta)
#print(sum((Y-f(sita))^2))
}


alpha=sita[1]
beta=sita[2]
Qinf=sita[3]

pred=Qinf/(1+alpha*exp(-beta*(t-t0)))

plot(t-t0,pred,xlim=c(min(t-t0),max(t-t0)),ylim=c(min(Q),max(Q)),xlab="t-t0",ylab="growth curve",col=2)
par(new=T)
plot(t-t0,Q,xlim=c(min(t-t0),max(t-t0)),ylim=c(min(Q),max(Q)),xlab="t-t0",ylab="growth curve",col=3)

```

```r

#p.202 

data(EuStockMarkets)
Q0=1670
Q=cumsum(as.numeric(t(EuStockMarkets[,2])))
Qdash=c(Q0,Q[1:(length(Q)-1)])
qt=Q-Qdash
qtdash=qt[1:(length(qt)-1)]
qt=qt[2:length(qt)]
t=(as.numeric(EuStockMarkets[,1]))/1000
t=t[2:(length(t))]
X=matrix(as.numeric(EuStockMarkets[,3:4]),nrow=nrow(EuStockMarkets))
X=apply(X,2,cumsum)
Xdash=X[1:(nrow(X)-1),]
X=X[2:nrow(X),]
t0=0


f=function(x){
b0=x[1];lam=x[2]
b1=x[-c(1:2)]
value=qtdash*lam+(b0+X%*%b1)/(1+alpha*exp(-beta*t))-(b0+Xdash%*%b1)/(1+alpha*exp(-beta*tdash))
return(value)
}


Y=qt
h=0.001
lam=1
ite=10^3
sita=c(5,1,5,5)
eta=10^(-2)

for(l in 1:ite){
sita_pre=sita
Z=array(0,dim=c(length(Y),length(sita)))
vec=sita
for(j in 1:ncol(Z)){
vec_sub=vec;vec_sub[j]=vec_sub[j]+h
Z_sub=(f(vec_sub)-f(vec))/h
Z[,j]=Z_sub
}
bb=solve(t(Z)%*%Z+diag(lam,ncol(Z)))%*%t(Z)%*%(Y-f(sita))
sita=c(sita+eta*bb)
print(sum((Y-f(sita))^2))
}

plot(t,f(sita),xlim=c(min(t),max(t)),ylim=c(min(Y),max(Y)),xlab="t",ylab="qt",col=2)
par(new=T)
plot(t,Y,xlim=c(min(t),max(t)),ylim=c(min(Y),max(Y)),xlab="t",ylab="qt",col=3)


```

```r

#p.202 

data(EuStockMarkets)
Q0=1670
Q=cumsum(as.numeric(t(EuStockMarkets[,2])))
Qdash=c(Q0,Q[1:(length(Q)-1)])
qt=Q-Qdash
qtdash=qt[1:(length(qt)-1)]
qt=qt[2:length(qt)]
t=(as.numeric(EuStockMarkets[,1]))/1000
tdash=t[1:(length(t)-1)]
t=t[2:(length(t))]
X=matrix(as.numeric(EuStockMarkets[,3:4]),nrow=nrow(EuStockMarkets))
X=apply(X,2,cumsum)
Xdash=X[1:(nrow(X)-1),]
X=X[2:nrow(X),]
t0=0


f=function(x){
b0=x[1];lam=x[2];alpha=x[3];beta=x[4]
b1=x[-c(1:4)]
value=qtdash*lam+(b0+X%*%b1)/(1+alpha*exp(-beta*t))-(b0+Xdash%*%b1)/(1+alpha*exp(-beta*tdash))
return(value)
}


Y=qt
h=0.001
lam=1
ite=4000
sita=c(5,1,5,5,5,5)
eta=10^(-2)

for(l in 1:ite){
sita_pre=sita
Z=array(0,dim=c(length(Y),length(sita)))
vec=sita
for(j in 1:ncol(Z)){
vec_sub=vec;vec_sub[j]=vec_sub[j]+h
Z_sub=(f(vec_sub)-f(vec))/h
Z[,j]=Z_sub
}
bb=solve(t(Z)%*%Z+diag(lam,ncol(Z)))%*%t(Z)%*%(Y-f(sita))
sita=c(sita+eta*bb)
sig=sum((Y-f(sita))^2)/(length(Y))
loglik=-length(Y)*log(2*pi*sig)/2-sum((Y-f(sita))^2)/(2*sig)
print(loglik)
#print(sum((Y-f(sita))^2))
}

plot(t,f(sita),xlim=c(min(t),max(t)),ylim=c(min(Y),max(Y)),xlab="t",ylab="qt",col=2)
par(new=T)
plot(t,Y,xlim=c(min(t),max(t)),ylim=c(min(Y),max(Y)),xlab="t",ylab="qt",col=3)


```