# CHAPTER 7: R COMMANDS # 08/01/2007 ################################################################################ # R PROGRAM TO ESTIMATE THE COEFFICIENTS OF AN AR(p) MODEL p=5 x=matrix(scan("Eurostoxx50.txt"),ncol=5,byrow=T)[,4] T=length(x) W=10000 mu=mean(x) sig2=var(x) lambdareal=2*runif(p)-1 preal=p pcomp=0 lambdacomp=0 llog=function(pr,pc,lr,lc,compsi=T,pepsi=0){ # LIKELIHOOD REPRESENTATION if (compsi) { Psi=matrix(0,ncol=p,nrow=p+1) Psi[1,]=1 if (pr>0) { Psi[2,1]=-lr[1] if (pr>1) { for (i in 2:pr) Psi[2:(i+1),i]=Psi[2:(i+1),i-1]-lr[i]*Psi[1:i,i-1] } } if (pc>0) { if (pr>0) { Psi[2,pr+2]=-2*lc[1]+Psi[2,pr] Psi[3:(pr+3),pr+2]=(lc[1]^2+lc[2]^2)*Psi[1:(pr+1),pr]-2*lc[1]*Psi[2:(pr+2),pr]+Psi[3:(pr+3),pr] } else { Psi[2,2]=-2*lc[1]; Psi[3,2]=(lc[1]^2+lc[2]^2); } if (pc>2) { for (i in seq(4,pc,2)) { pri=pr+i prim=pri-2 Psi[2,pri]=-2*lc[i-1]+Psi[2,prim] Psi[3:(pri+1),pri]=(lc[i-1]^2+lc[i]^2)*Psi[1:(pri-1),prim]-2*lc[i-1]*Psi[2:pri,prim]+Psi[3:(pri+1),prim] } } } Psi=Psi[1:(p+1),p] } else { Psi=pepsi } # LOGLIKELIHOOD x=x-mu loglike=-x[1]^2 for (i in 2:p) loglike=loglike-(t(Psi[1:i])%*%x[i:1])^2 for (i in (p+1):T) loglike=loglike-(t(Psi)%*%x[i:(i-p)])^2 loglike=(loglike/sig2-T*log(sig2))/2 x=x+mu list(ll=loglike,ps=Psi) } llo=llog(pr=preal,pc=pcomp,lr=lambdareal,lc=lambdacomp) # MCMC REVERSIBLE SCHEME psis=matrix(0,ncol=p,nrow=W) mus=rep(0,W) sigs=rep(0,W) ncomp=rep(0,W) llik=rep(0,W) indacpt=rep(0,5) for (m in 1:W){ if (runif(1)<.1{ ind=sample(1:p,1) if (ind<=pcomp){ ind=ind-(ind%%2==0) ppropreal=preal ppropcomp=pcomp lambpropreal=lambdareal lambpropcomp=lambdacomp lambpropcomp[ind]=lambdacomp[ind]+.05*rnorm(1) lambpropcomp[ind+1]=lambdacomp[ind+1]+.05*rnorm(1) }else{ ppropreal=preal ppropcomp=pcomp lambpropreal=lambdareal lambpropcomp=lambdacomp lambpropreal[ind-pcomp]=lambdareal[ind-pcomp]+.05*rnorm(1) } lloprop=llog(pr=ppropreal,pc=ppropcomp,lr=lambpropreal,lc=lambpropcomp) if (log(runif(1))1)&&(pcomp>0)){ if (runif(1)<.1){ ppropcomp=pcomp-2 ppropreal=preal+2 ind=sample(1:pcomp,1) ind=ind-(ind%%2==0) if (ppropcomp>0){ lambpropcomp=lambdacomp[((1:pcomp)!=ind)&((1:pcomp)!=(ind+1))] }else{ lambpropcomp=0 } lambpropreal=c(lambdareal,lambdacomp[ind]+0.05*rnorm(2)) coef=9*(1+(preal<2)) * (pi/4) }else{ ppropreal=preal-2 ppropcomp=pcomp+2 ind=sample(1:preal,2) #real roots to remove if (ppropreal>0){ lambpropreal=lambdareal[((1:preal)!=ind[1])&((1:preal)!=ind[2])] }else{ lambpropreal=0 } lambpropcomp=c(lambdacomp,mean(lambdareal[ind])+.05*rnorm(1),.05*rnorm(1)) coef=(4/pi) * (1+(ppropcomplloprop\$ll-llo\$ll){ indacpt[3]=1+indacpt[3] mu=muold }else{ llo=lloprop } sig2old=sig2 sig2=exp(rnorm(1,mean=log(sig2),sd=sqrt(2)*sig2)) lloprop=llog(pr=preal,pc=pcomp,lr=lambdareal,lc=lambdacomp,compsi=F,pepsi=llo\$ps) if (log(runif(1))>lloprop\$ll-llo\$ll){ indacpt[4]=1+indacpt[4] sig2=sig2old }else{ llo=lloprop } } psis[m,]=llo\$ps[2:(p+1)] mus[m]=mu sigs[m]=sig2 llik[m]=llo\$ll ncomp[m]=pcomp } # GRAPHICAL REPRESENTATION FOR THE AR(p) MODEL par(mfrow=c(3,3),mar=c(4,4,2,1)) if (p>3) { hist(ncomp,main="",xlab="p",ylab="",col="gold4",breaks=seq(-1,p+1,2));par(new=T);plot(ncomp,axes=F,cex=.3,xlab="",ylab="") axis(side=4) } if (p<4) { plot(ncomp,cex=.3,xlab="Iterations",ylab="Complex roots") } plot(mus,type="l",col="steelblue4",xlab="Iterations",ylab=expression(mu)) plot(500:W,sigs[500:W],type="l",col="steelblue4",xlab="Iterations",ylab=expression(sigma^2)) for (i in 1:min(3)) plot(500:W,psis[500:W,i],type="l",col="steelblue4",xlab="Iterations",ylab=expression(psi)) plot(llik,type="l",col="sienna4",xlab="Iterations",ylab="log-likelihood") pst=matrix(1,ncol=(p+1),nrow=W);pst[,2:(p+1)]=psis lame=apply(pst,1,polyroot) plot((1/lame)[Mod(lame)>1],col="gold",cex=.3,xlab=expression(Re(lambda)),ylab=expression(Im(lambda))) lines(seq(-1,1,.01),sqrt(1-seq(-1,1,.01)^2),col="sienna",lty=2,lwd=2) lines(seq(-1,1,.01),-sqrt(1-seq(-1,1,.01)^2),col="sienna",lty=2,lwd=2) pses=apply(psis,2,mean) mimine=mean(mus) predo=mimine-pses[1]*(x[p:(T-1)]-mimine) for (i in 2:p) predo=predo-pses[i]*(x[(p-i+1):(T-i)]-mimine) plot(x[(p+1):T],type="l",col="steelblue4",xlab="t",ylab="x") lines(predo,lty=2,col="sienna4",lwd=1.8) ################################################################################ # R PROGRAM TO ESTIMATE THE COEFFICIENTS OF AN AR(q) MODEL p=9 q=p T=550 x=matrix(scan("Dyndata"),ncol=5,byrow=T)[1:(T+50),5] # observations x=x[1:T] W=10000 mu=mean(x) varef=var(x) sig2=varef lambdareal=2*runif(p)-1 preal=p pcomp=0 lambdacomp=0 eps=rnorm(p,sd=sqrt(sig2)) #Just about anythin' llog=function(pr,pc,lr,lc,compsi=T,pepsi=0,eps) { # LIKELIHOOD REPRESENTATION if (compsi) { Psi=matrix(0,ncol=p,nrow=p+1) Psi[1,]=1 if (pr>0) { Psi[2,1]=-lr[1] if (pr>1) { for (i in 2:pr) Psi[2:(i+1),i]=Psi[2:(i+1),i-1]-lr[i]*Psi[1:i,i-1] } } if (pc>0) { if (pr>0) { Psi[2,pr+2]=-2*lc[1]+Psi[2,pr] Psi[3:(pr+3),pr+2]=(lc[1]^2+lc[2]^2)*Psi[1:(pr+1),pr]-2*lc[1]*Psi[2:(pr+2),pr]+Psi[3:(pr+3),pr] } else { Psi[2,2]=-2*lc[1]; Psi[3,2]=(lc[1]^2+lc[2]^2); } if (pc>2) { for (i in seq(4,pc,2)) { pri=pr+i prim=pri-2 Psi[2,pri]=-2*lc[i-1]+Psi[2,prim] Psi[3:(pri+1),pri]=(lc[i-1]^2+lc[i]^2)*Psi[1:(pri-1),prim]-2*lc[i-1]*Psi[2:pri,prim]+Psi[3:(pri+1),prim] } } } Psi=Psi[2:(p+1),p] } else { Psi=pepsi } x=x-mu heps=rep(0,T+p) heps[1:p]=eps for (i in 1:T) heps[p+i]=x[i]+sum(rev(Psi)*heps[i:(p+i-1)]) loglike=-((sum(heps^2)/sig2)+(T+p)*log(sig2))/2 x=x+mu list(ll=loglike,ps=Psi) } llo=llog(pr=preal,pc=pcomp,lr=lambdareal,lc=lambdacomp,eps=eps) # MCMC REVERSIBLE SCHEME psis=matrix(0,ncol=p,nrow=W) mus=rep(0,W) sigs=rep(0,W) ncomp=rep(0,W) llik=rep(0,W) indacpt=rep(0,10) epsrec=matrix(0,ncol=p,nrow=W) preds=matrix(0,ncol=p,nrow=W) for (m in 1:W){ ind=sample(1:p,1) if (ind<=pcomp){ ind=ind-(ind%%2==0) ppropreal=preal ppropcomp=pcomp lambpropreal=lambdareal lambpropcomp=lambdacomp lambpropcomp[ind]=lambdacomp[ind]+.05*rnorm(1) lambpropcomp[ind+1]=lambdacomp[ind+1]+.05*rnorm(1) }else{ ppropreal=preal ppropcomp=pcomp lambpropreal=lambdareal lambpropcomp=lambdacomp lambpropreal[ind-pcomp]=lambdareal[ind-pcomp]+.05*rnorm(1) } lloprop=llog(pr=ppropreal,pc=ppropcomp,lr=lambpropreal,lc=lambpropcomp,eps=eps) if (log(runif(1))1)&&(pcomp>0)){ if (runif(1)<.1){ ppropcomp=pcomp-2 ppropreal=preal+2 ind=sample(1:pcomp,1) ind=ind-(ind%%2==0) if (ppropcomp>0){ lambpropcomp=lambdacomp[((1:pcomp)!=ind)&((1:pcomp)!=(ind+1))] }else{ lambpropcomp=0 } lambpropreal=c(lambdareal,lambdacomp[ind]+0.05*rnorm(2)) coef=9*(1+(preal<2)) * (pi/4) }else{ ppropreal=preal-2 ppropcomp=pcomp+2 ind=sample(1:preal,2) #real roots to remove if (ppropreal>0){ lambpropreal=lambdareal[((1:preal)!=ind[1])&((1:preal)!=ind[2])] }else{ lambpropreal=0 } lambpropcomp=c(lambdacomp,mean(lambdareal[ind])+.05*rnorm(1),.05*rnorm(1)) coef=(4/pi) * (1+(ppropcomplloprop\$ll-llo\$ll){ indacpt[3]=1+indacpt[3] mu=muold }else{ llo=lloprop } if (runif(1)<.3){ Psi=llo\$ps heps=rep(0,2*p) heps[1:p]=eps for (j in (p+1):(2*p)) heps[j]=x[j]-sum(Psi*heps[(j-p):(j-1)]) varheps=(2*p)*var(heps) sig2old=sig2 sig2=varheps/rgamma(1,2*p) difdens=dgamma(sig2/varheps,2*p,log=T)-dgamma(sig2old/varheps,2*p,log=T)+log(sig2)-log(sig2old) lloprop=llog(pr=preal,pc=pcomp,lr=lambdareal,lc=lambdacomp,compsi=F,pepsi=llo\$ps,eps=eps) if (log(runif(1))>lloprop\$ll-llo\$ll-difdens){ sig2=sig2old }else{ indacpt[4]=3+indacpt[4] llo=lloprop } } if (runif(1)<.3){ sig2old=sig2 sig2=exp(rnorm(1,mean=log(sig2),sd=sqrt(.1*varef))) lloprop=llog(pr=preal,pc=pcomp,lr=lambdareal,lc=lambdacomp,compsi=F,pepsi=llo\$ps,eps=eps) difdens=log(sig2)-log(sig2old) if (log(runif(1))>lloprop\$ll-llo\$ll-difdens){ sig2=sig2old }else{ llo=lloprop indacpt[5]=3+indacpt[5] } } if (runif(1)<.3){ sig2old=sig2 thismean=varef/(1+sum(Psi^2)) sig2=exp(rnorm(1,mean=thismean,sd=sqrt(varef))) difdens=dnorm(log(sig2),mean=thismean,sd=sqrt(.5*varef),log=T)-dnorm(log(sig2old),mean=thismean,sd=sqrt(.5*varef),log=T) lloprop=llog(pr=preal,pc=pcomp,lr=lambdareal,lc=lambdacomp,compsi=F,pepsi=llo\$ps,eps=eps) if (log(runif(1))>lloprop\$ll-llo\$ll-difdens){ sig2=sig2old }else{ llo=lloprop indacpt[6]=3+indacpt[6] } } Psi=llo\$ps if (runif(1)<.5){ heps=rep(0,2*p+1) keps=rep(0,p) for (i in 1:q){ x = x-mu heps[1:p]=eps for (j in (p+1):(2*p+1)) heps[j]=x[j]+sum(rev(Psi)*heps[(j-p):(j-1)]) heps[i]=0 for (j in 1:(q-i+1)) keps[j]=x[j]+sum(rev(Psi)*heps[j:(j+p-1)]) x = x+mu epsvar = 1/sum(c(1,Psi[i:q]^2)) epsmean = sum(Psi[i:q]*keps[1:(q-i+1)])*epsvar epsmean = epsmean/epsvar epsvar = sig2*epsvar propeps = rnorm(1,mean=epsmean,sd=sqrt(epsvar)) epspr=eps epspr[i]=propeps lloprop=llog(pr=preal,pc=pcomp,lr=lambdareal,lc=lambdacomp,compsi=F,pepsi=Psi,eps=epspr) propsal1=dnorm(propeps,mean=epsmean,sd=sqrt(epsvar),log=T) x = x-mu heps[i]=propeps for (j in (p+1):(2*p+1)) heps[j]=x[j]+sum(rev(Psi)*heps[(j-p):(j-1)]) heps[i]=0 for (j in 1:(q-i+1)) keps[j]=x[j]+sum(rev(Psi)*heps[j:(j+p-1)]) x = x+mu epsvar = 1/sum(c(1,Psi[i:q]^2)) epsmean = sum(Psi[i:q]*keps[1:(q-i+1)]) epsmean = epsmean*epsvar epsvar = sig2*epsvar propsal0=dnorm(eps[i],mean=epsmean,sd=sqrt(epsvar),log=T) if (log(runif(1))1000) y=y[seq(1,t[1],length=1000),] if (t[2]>1000) y=y[,seq(1,t[2],length=1000)] } else y=y[seq(1,length(x),length=1000)] y } start=1000 par(mfrow=c(3,3),mar=c(4,4,2,1)) if (p>3) { hist(ncomp,main="",xlab="p",ylab="",col="gold4",breaks=seq(-1,p+1,2));par(new=T); plot(sube(1:W),sube(ncomp),axes=F,cex=.3,xlab="",ylab="") axis(side=4) } if (p<4) { plot(sube(1:W),sube(ncomp),cex=.3,xlab="Iterations",ylab="Complex roots") } plot(sube(start:W),sube(mus[start:W]),type="l",col="steelblue4",xlab="Iterations",ylab=expression(mu)) plot(sube(start:W),sube(sigs[start:W]),type="l",col="steelblue4",xlab="Iterations",ylab=expression(sigma^2)) for (i in 1:min(3)) plot(sube(start:W),sube(psis[start:W,i]),type="l",col="steelblue4",xlab="Iterations",ylab=expression(psi)) plot(sube(start:W),sube(llik[start:W]),type="l",col="sienna4",xlab="Iterations",ylab="log-likelihood") pst=matrix(1,ncol=(p+1),nrow=W);pst[,2:(p+1)]=-psis lame=apply(pst,1,polyroot) plot(sube((1/lame)[Mod(lame)>1]),col="gold",cex=.3,xlab=expression(Re(lambda)),ylab=expression(Im(lambda))) lines(seq(-1,1,.01),sqrt(1-seq(-1,1,.01)^2),col="sienna",lty=2,lwd=2) lines(seq(-1,1,.01),-sqrt(1-seq(-1,1,.01)^2),col="sienna",lty=2,lwd=2) plot(sube(start:W),sube(epsrec[start:W,1]),ylim=range(epsrec[start:W,]),type="l",ylab=expression(epsilon),col="steelblue") for (i in 2:p) lines(sube(start:W),sube(epsrec[start:W,i]),col="steelblue")