***************************************************************** * * * Ref: Pai and Ravishanker, JTSA 1998 * ***************************************************************** * modified gibbs sampling for arfima(p,d,q) models with mean * * uses exact likelihood with np+nq latent variables * ***************************************************************** ***************************************************************** * The model fitted is * ***************************************************************** * PHI(b)*(1-b)**d*(y(t)-mu) = THETA(b)*e(t) * * where * * PHI(b) = (1 - phi(1)*b - ... - phi(p)*b**p), * * THETA(b) = (1 - theta(1)*b - ... - theta(q)*b**q), * * y(t) is the given time series, * * e(t) is the white noise series which is normal(0,vari) * * mu is the mean of the series, * * d is the fractional differencing parameter * ***************************************************************** ***************************************************************** * Metropolis-within-Gibbs * ***************************************************************** * This program uses metropolis-within-gibbs * * to generate samples * * from the posterior distributions of * * PHI: vector of dim. np of ar parameters * * THETA: vector of dim. nq of ma parameters * * d: fractional differencing parameter * * mu: scalar mean parameter * * yo: vector of dim. np of latent variables * * unobserved historical data * * a0: vector of dim. nq of latent variables * * unobserved fractional noise * * sigma: scalar white noise standard deviation * ***************************************************************** ***************************************************************** * Stationary and invertible regions * ***************************************************************** * PHI in Cp, THETA in Cq and d in (-.5,.5). * ***************************************************************** ***************************************************************** * Gibbs blocking * ***************************************************************** * (PHI,THETA,d): Metropolis with truncated Normal probing* * (y0,a0): Metropolis with Normal probing dist. * * mu: Metropolis with Normal probing dist. * * sigma: Inverse Gamma * ***************************************************************** ***************************************************************** * Description of parameter inputs * ***************************************************************** * np: number of ar parameters. * * nq: number of ma parameters. * * ixd: ixd=0 if short memory (arma) models * * ixd=1 if long memory (arfima) models * * npq: np+nq * * npqd: npq+ixd * * n: number of data used for the model fitting * * nnq: n+nq * * iter1: number of iterations in the gibbs sampler * * nrep1: number of replications in the gibbs sampler * * nloop: number of trajectories in the metropolis chains * * nw: # of terms used in hosking's (1981, eqn 4.1)) * * to compute the acf of arfima(p,d,q) process * * nmx: size of work array, >= n * * npqmx: size of work array, >= np+nq+ixd * * iterx: size of work array, >= iter1 * * nrepx: size of work array, >= nrep1 * * nloopx: size of work array, >= nloop * * nmx2: size of work array, * * nmx2 >= (n+nq+1)**2/4 if n+nq is odd * * nmx2 >= (n+nq+2)*(n+nq)/4 if n+nq is even * ***************************************************************** * e.g. n=1000,nq=10 => nmx=255530 * e.g. n=400 ,nq=10 => nmx=42230 ***************************************************************** ***************************************************************** * variables in the program * ***************************************************************** * zdata(n): total observations * * data(n) : mean subtracted data * * rsig(npqd*(npqd+1)/2): VAR(PHI,THETA,d) & its Cholesky * * for Metropolis probing * * omgch(npq*(npq+1)/2): VAR(y0,a0) & its Cholesky Deco. * * for Metropolis probing * ***************************************************************** * output from Gibbs sampler * ***************************************************************** * mphi1(iter1,nrep1,np): ar parameter vector * * mthe1(iter1,nrep1,nq): ma parameter vector * * mean1(iter1,nrep1): mean parameter * * sigma1(iter1,nrep1): white noise std. dev. * * y01(iter1,nrep1,1-np:0): latent variables * * (unobserved data) * * eps01(iter1,nrep1,1-nq:0): latent variables * * (unobserved fractional noise) * * rd1(iter1,nrep1): fractional difference parameter * ***************************************************************** ***************************************************************** * routines called and purpose: * ***************************************************************** * * * lm: main program, can be called from Splus * * * * armaco: covariance or correlation of an arma model * * includes dterm and lineq * * chod: compute cholesky decomposition of a matrix * * fcal: compute log of the posterior density * * to which the complete conditional distributions * * are proportional * * dmvnor: random sample from a multivariate normal dist. * * g1: compute the sum of squares * * in the function g1 (see pai and ravishanker) * * g2: compute the sum of squares & determinant * * in the function g2 (see pai and ravishanker) * * gammad: random sample from a Gamma dist. * * includes gamma1,gamma2 and expon * * gama: gamma function * * gc: compute the determinant and * * inverse of the covariance matrix of dimension * * (n+nq,n+nq) from arfima(0,d,0) process * * gcova: computes the acf of an arfima(0,d,0) process * * gctran: supressed a double symmetric matrix * * gensig: generate sample for white noise std. dev. * * gsy: computes the acf of an arfima(p,d,q) proces * * and eta components * * musig: compute the runeps0 vector (of dim. 1-nq:n) * * from arfima(0,d,0) process * * perturb: generate starting values for parallel gibbs * * chains by perturbing the initial value * * rand: random sample from an uniform dist. * * troot: test stationarity or invertibility * * upeps: generate sample for (eps0,y0), latent variables * * update: generate sample for (PHI,THETA,d) * * updatemu: generate sample for mean mu * * unormd: random sample from an univariate normal dist. * * xomega: compute covariance matrix of (y0,a0) * * (similar to newbold, biometrika,1974) * * and its Cholesky decomposition * * xsminv: inverts a symmetric matrix * ***************************************************************** ***************************************************************** * SUBROUTINE lm1 * ***************************************************************** subroutine lm1 (datain,nin,para,npara,ap1,nap1,ap2,nap2, 6 mphi1,mthe1,rd1,mean1,sigma1,y01,eps01) ***************************************************************** parameter (npqmx=10,nmx=400) parameter (nmx2=42230) parameter (iterx=9000) parameter (nrepx=10) parameter (nloopx=20) parameter (nw=500) implicit real*8 (a-h,o-z) real*8 datain(nin) real*8 ap1(nap1),ap2(nap2) real*8 sigma1(iterx,nrepx) real*8 zdata(nmx),data(nmx) real*8 runphi(0:npqmx),runthe(0:npqmx) real*8 mphi1(iterx,nrepx,npqmx) real*8 mthe1(iterx,nrepx,npqmx) real*8 mean1(iterx,nrepx) real*8 rd1(iterx,nrepx) real*8 y01(iterx,nrepx,1-npqmx:0) real*8 eps01(iterx,nrepx,1-npqmx:0) real*8 runy0(1-npqmx:0),runeps0(1-npqmx:0) real*8 runeps(1-npqmx:nmx) real*8 rsig(npqmx*(npqmx+1)/2) real*8 omgch(npqmx*(npqmx+1)/2) real*8 c1(nmx2) real*8 cova(0:nw),covy(0:npqmx-1) real*8 s(1-npqmx-nmx:npqmx-1) integer para(npara),ixd common/nconst/np,nq,ixd,npq,npqd,n,nnq,iter1,nrep1,nloop ***************************************************************** * Set up parameter for this subroutine * ***************************************************************** np= para(1) nq= para(2) ixd= para(3) iter1= para(4) nrep1= para(5) nloop= para(6) c npq= np+nq npqd= npq+ixd n= nin nnq= n+nq c stdmu=ap1(npqd+2) runphi(0)=-1.0d0 runthe(0)=-1.0d0 ***************************************************************** * input data in zdata * ***************************************************************** do 10 i=1,nin zdata(i)=datain(i) 10 continue ***************************************************************** * input covariance matrix of (PHI,THETA,d) and * * convert rsig to its cholesky decomposition * ***************************************************************** do 40 i=1,nap2 40 rsig(i)=ap2(i) call chod (rsig,npqd) ***************************************************************** * form the covariance matrix of e*(latent variable) * * and then its cholesky decomposition in omgch * ***************************************************************** if (npq.ne.0) call xomega (omgch,ap1,nap1) ***************************************************************** * perturb the initial values to obtain * * starting values for nrep1 parallel gibbs chains. * ***************************************************************** call perturb (mphi1,mthe1,y01,eps01,sigma1,mean1,rd1 / ,rsig,omgch,ap1,nap1) time_now = wall_seconds() ***************************************************************** * begin nrep1 independent Gibbs chains * ***************************************************************** do 1120 jj=1,nrep1 write (*,99) 'Replication #:', $ jj,wall_seconds()-time_now 99 format(a,2x,i5,2x,f15.3) ***************************************************************** * set starting value of each chain * ***************************************************************** if (np.ne.0) then do 1126 i=1,np 1126 runphi(i)=mphi1(1,jj,i) endif c if (nq.ne.0) then do 1128 i=1,nq 1128 runthe(i)=mthe1(1,jj,i) endif c vari=sigma1(1,jj)**2 runmn=mean1(1,jj) if (ixd.eq.1) rd=rd1(1,jj) ***************************************************************** * subtracts mean from zdata to form data * ***************************************************************** do 1121 i=1,n 1121 data(i)=zdata(i)-runmn c if (nq.ne.0) then do 1122 i=1,nq 1122 runeps0(i-nq)=eps01(1,jj,i-nq) endif c if (np.ne.0) then do 1124 i=1,np 1124 runy0(i-np)=y01(1,jj,i-np)-runmn endif ***************************************************************** call musig (runphi,runthe,runy0,runeps0,data,runeps) ***************************************************************** call gc (rd,de1,c1) call g1 (runeps,c1,a1) if (np.ne.0) then if (ixd.eq.1) call gcova (cova,rd) call gsy (cova,runphi,runthe,s,covy) call g2 (c1,runeps,runy0,covy,s,rd,a2,de2) endif ***************************************************************** * run MCMC with iter1-1 iteration * ***************************************************************** do 1110 ii=1,iter1-1 write (*,99)'MC Iteration: ', ii, $ wall_seconds()-time_now call fcal (f,a1,a2,de1,de2,vari) c write(6,*)'f,a1,a2,de1,de2',f,a1,a2,de1,de2 ***************************************************************** * update latent variables * ***************************************************************** if (npq.ne.0) then call upeps (f,c1,de1,runphi,runthe,runy0,runeps0, / omgch,data,runeps,vari,rd,cova,covy,s) c if (nq.ne.0) then do 1134 k=1,nq 1134 eps01(ii+1,jj,k-nq)=runeps0(k-nq) endif c if (np.ne.0) then do 1132 k=1,np 1132 y01(ii+1,jj,k-np)=runy0(k-np)+runmn endif endif c c c write(6,*)'after upeps f a1 a2 de1 de2',f,a1,a2,de1,de2 ***************************************************************** * update (PHI,THETA,d) * ***************************************************************** call update (f,c1,de1,rd,runphi,runthe,runy0,runeps0, / vari,a1,a2,rsig,runeps,data,cova,covy,s) c if (np.ne.0) then do 1140 i=1,np 1140 mphi1(ii+1,jj,i)=runphi(i) endif c if (nq.ne.0) then do 1150 i=1,nq 1150 mthe1(ii+1,jj,i)=runthe(i) endif c if (ixd.eq.1) rd1(ii+1,jj)=rd c c write(6,*)'after update f,a1,a2,de1,de2',f,a1,a2,de1,de2 ***************************************************************** * update mu * ***************************************************************** call updatemu (f,c1,de1,runmn,runphi,runthe,stdmu, / data,runy0,runeps0,runeps,vari,rd,cova,covy,s) mean1(ii+1,jj)=runmn c c write(6,*)'after updatemu f,a1 a2 de1,de2',f,a1,a2,de1,de2 ***************************************************************** * update sigma * ***************************************************************** call gensig (a1,a2,vari) sigma1(ii+1,jj)=dsqrt(vari) c ***************************************************************** 1110 continue 1120 continue call prx0 (mphi1,mthe1,y01,eps01,mean1,sigma1,rd1) return end ********************************************************************** * SUBROUTINE chod * ********************************************************************** subroutine chod (x,n) parameter(npqmx=10) implicit real*8 (a-h,o-z) real*8 x(npqmx*(npqmx+1)/2) ********************************************************************** c* * c* purposes: function xsminv inverts a given (n*n)-dimensioned * c* symmetric positive definite matrix, using the * c* cholesky's upper triangular factorization method. * c* * c* usage : det=xsminv(y,n) * c* * c* * c* parameters on entry on return * c* * c* x : corvariance its decomposition * c* * c* n : number of rows unchanged * c* * c* remark : error glag xsminv is set to the following values: * c* * c* xsminv > 0.0 : no error - matrix's determination * c* xsminv = 0.0 : singular matrix * c* xsminv = -1.0 : wrong parameter n on entry * c* xsminv = -5.0 : matrix not definite positive * c* * ********************************************************************** c det=1.0 if (n.le.0) then xsminv=-1.0 return endif c c compute upper-triangular factorization matrix t, c such as x=t'*t, by the cholesky's square root method. c do 1000 i=1,n max=(i*i+i)/2 ijk=max do 1000 j=i,n y=0.0 if (i.gt.1) then do 100 k=1,i-1 l=max-k m=ijk-k 100 y=y+x(l)*x(m) endif y=x(ijk)-y if (j.eq.i) then if (y.lt.0.0) then xsminv=-5.0 return elseif (y.eq.0.0) then xsminv=0.0 return endif z=dsqrt(y) x(max)=z det=det*x(max) z=1.0/z else x(ijk)=y*z endif 1000 ijk=ijk+j c return end ********************************************************************** * * * SUBROUTINE dmvnor * * * ********************************************************************** subroutine dmvnor (chod1,rn,n) ***************************************************************** * generate random sample rn(n) from multivariate Normal * * with mean 0 and chod1 as Cholasky decomposition * * of its variance * ***************************************************************** implicit real*8 (a-h, o-z) parameter (npqmx=10) dimension chod1(npqmx*(npqmx+1)/2),rn(npqmx) c do 5 i=1,n 5 rn(i)=unormd() c do 10 i=1,n a1=0.0d0 do 20 j=i,n 20 a1=a1+rn(j)*chod1(j*(j-1)/2+i) 10 rn(i)=a1 c return end ********************************************************************** * * * SUBROUTINE fcal * * * ********************************************************************** subroutine fcal (f,a1,a2,de1,de2,vari) ********************************************************************** * use: compute the log of the posterior density * * parameters of routine: * * f *output* scalar.log of the posterior density * * a1 *input * scalar.sum of squares term in the g1 function * * a2 *input * scalar. sum of squares term in the g2 function * * de1 *input * scalar. determinant in the g1 function * * de2 *input * scalar. determinant in the g2 function * * vari *input * scalar. contains white noise variance * ********************************************************************** parameter (npqmx=10,nmx=400) implicit real*8 (a-h,o-z) common/nconst/np,nq,ixd,npq,npqd,n,nnq,iter1,nrep1,nloop c c dpi=3.14159265359d0 c if (np.eq.0) then f=-0.5d0*(dble(npq+n+2)*dlog(2*vari)+de1+a1/vari) else f=-0.5d0*(dble(npq+n+2)*dlog(2*vari)+de1+de2+(a1+a2)/vari) endif c return end ********************************************************************** * * * SUBROUTINE g1 * * * ********************************************************************** subroutine g1 (runeps,c1,a1) ********************************************************************** * use: compute sum of squares term in g1 function * * parameters of routine: * * runeps *input * array of length (1-nq:n). * * c1 *input * array of length nnq*(nnq+1)/2. contains elements * * of the inverse of sigma22 * * a1 *output* scalar. computed value of sum of squares term * * in the g1 function * ********************************************************************** parameter (nmx=400,npqmx=10) parameter (nmx2=42230) implicit real*8 (a-h,o-z) real*8 runeps(1-npqmx:nmx) real*8 c1(nmx2) common/nconst/np,nq,ixd,npq,npqd,n,nnq,iter1,nrep1,nloop c c nnq=n+nq a1=0.0d0 c if (ixd.eq.0) then do 340 i=1,nnq a1=a1+runeps(i-nq)**2 340 continue else do 345 i=1,nnq do 345 j=1,nnq call gctran (i,j,ij,nnq) 345 a1=a1+c1(ij)*runeps(i-nq)*runeps(j-nq) endif return end ********************************************************************** * * * SUBROUTINE g2 * * * ********************************************************************** subroutine g2 (c1,runeps,runy0,covy,s,rd,a2,de2) ********************************************************************** * use: compute determinant and sum of squares in the g2 function * * i.e. determinant and inverse of sigma22 * * parameters of routine: * * c1 *input * array of length nnq*(nnq+1)/2 contains elements * * of inverse of sigma22 * * runeps *input * array of length 1-nq:0. contains current eps0 * * runy0 *input * array of length 1-np:0. contains current y0 * * s *input * array of length 1-np-n:nq-1. eta(k) elements * * rd *input * scalar. contains current value of fractional * * difference parameter d * * a2 *output* scalar. sum of squares in g2 function * * de2 *output* scalar. determinant in g2 function * * other routines used: xsminv * ********************************************************************** parameter (nmx=400,npqmx=10) parameter (nmx2=42230) parameter (nw=500) implicit real*8 (a-h,o-z) real*8 runeps(1-npqmx:nmx),runy0(1-npqmx:0) real*8 c1(nmx2) real*8 covy(0:npqmx-1) real*8 s(1-npqmx-nmx:npqmx-1) real*8 sig112(npqmx,npqmx),vmu(npqmx) real*8 sigt(npqmx,nmx+npqmx) common/nconst/np,nq,ixd,npq,npqd,n,nnq,iter1,nrep1,nloop ********************************************************************** * short memory ar(p) models * ********************************************************************** if (ixd.eq.0 .and. nq.eq.0) then do 400 i=1,np do 400 j=i,np sig112(i,j)=covy(j-i) 400 if (i.ne.j) sig112(j,i)=sig112(i,j) c det=xsminv(sig112,np) c if (det.le.0.0) then a2=1.0d10 de2=1.0d10 call dblepr (" 6 sigma112 is not positive definite det=",50,det,1) ********12345678901234567890123456789012345678901234567890" return else de2=dlog(det) endif c a2=0.0d0 do 420 i=1,np do 420 j=1,np 420 a2=a2+sig112(i,j)*runy0(i-np)*runy0(j-np) return endif ********************************************************************** * short memory arma(p,q) models * ********************************************************************** if (ixd.eq.0 .and. nq.ne.0) then do 500 i=1,np do 500 j=i,np sig112(i,j)=covy(j-i) do 501 i1=1,nq 501 sig112(i,j)=sig112(i,j)-s(j-np-i1+nq)*s(i-np-i1+nq) 500 if (i.ne.j) sig112(j,i)=sig112(i,j) c do 510 i=1,np vmu(i)=runy0(i-np) do 510 i1=1,nq 510 vmu(i)=vmu(i)-s(i-np-i1+nq)*runeps(i1-nq) c det=xsminv(sig112,np) c if (det.le.0.0) then a2=1.0d10 de2=1.0d10 call dblepr (" 6 sigma112 is not pd,!!!!!!!!!!!!!!!!!!!!!!!!!!!det=",50,det,1) return else de2=dlog(det) endif c c a2=0.0d0 do 520 i=1,np do 520 j=1,np 520 a2=a2+vmu(i)*sig112(i,j)*vmu(j) c return endif ********************************************************************** * long memory arfima(p,d,q) models * ********************************************************************** if (ixd.eq.1) then c nnq=n+nq c do 100 i=1,np do 100 j=1,nnq sigt(i,j)=0. do 110 k=1,nnq call gctran (k,j,kj,nnq) 110 sigt(i,j)=sigt(i,j)+s(i-np-k+nq)*c1(kj) 100 continue c do 600 i=1,np do 600 j=i,np sig112(i,j)=covy(j-i) do 601 i1=1,nnq 601 sig112(i,j)=sig112(i,j)-sigt(i,i1)*s(j-np-i1+nq) 600 if (i.ne.j) sig112(j,i)=sig112(i,j) c do 610 i=1,np vmu(i)=runy0(i-np) do 610 i1=1,nnq 610 vmu(i)=vmu(i)-sigt(i,i1)*runeps(i1-nq) c det=xsminv(sig112,np) c if (det.le.0.0) then a2=1.0d10 de2=1.0d10 call dblepr (" 6 sigma112 is not pd,!!!!!!!!!!!!!!!!!!!!!!!!!!!det=",50,det,1) return else de2=dlog(det) endif c a2=0.0d0 do 620 i=1,np do 620 j=1,np 620 a2=a2+vmu(i)*sig112(i,j)*vmu(j) endif c return end ********************************************************************** * * * FUNCTION gammad * * * ********************************************************************** function gammad (alpha, beta) ************************************************* * * * This function generates a random gamma * * number, calling either gamma1 (alpha > 1.0) * * or expon (alpha = 1.0) or gamma2 (alpha < * * 1.0). * * * * Density is proportional to * * x**(alpha - 1)*exp(-x/beta) * * * * Reference RIPLEY, B.D. (1987) * * "STOCHASTIC SIMULATION" * * * ************************************************* implicit real*8 (a-h, o-z) if(alpha.gt.1.0d0) then gammad = gamma1(alpha, beta) else if (alpha.lt.1.0d0) then gammad = gamma2(alpha, beta) else gammad = expon(beta) end if return end ************************************************* * * * FUNCTION gamma1(alpha, beta) * * * ************************************************* function gamma1(alpha, beta) implicit real*8 (a-h, o-Z) save aprev, c1, c2, c3, c4, c5 data aprev /0.0d0/ if(alpha .eq. aprev) goto 10 c1 = alpha - 1.0d0 aa = 1.0d0/c1 c2 = aa*(alpha - 1.0d0/(6.0d0*alpha)) c3 = 2.0d0*aa c4 = c3 + 2.0d0 if(alpha .gt. 2.5d0) c5 = 1.0d0/dsqrt(alpha) 10 u1 = rand() u2 = rand() if(alpha .le. 2.5d0) goto 20 u1 = u2 + c5*(1.0d0 - 1.86*u1) if(u1 .le. 0.0d0 .or. u1 .ge. 1.0d0) goto 10 20 w = c2*u2/u1 if(c3*u1 + w + 1.0d0/w .lt. c4) goto 30 if(c3*dlog(u1) - dlog(w) + w .ge. 1.0d0) goto 10 30 gamma1 = c1*w*beta aprev = alpha return end ************************************************* * * * FUNCTION gamma2(alpha, beta) * * * ************************************************* function gamma2(alpha, beta) implicit real*8 (a-h, o-z) e = dexp(1.0d0) b = (alpha + e)/e 10 p = b*rand() if(p .gt. 1.0d0) goto 20 x = p**(1.0d0/alpha) ux = rand() if(x .gt. -dlog(ux)) goto 10 gamma2 = x*beta return 20 x = -dlog( (b-p)/alpha ) u1 = rand() if(x**(alpha - 1.0d0) .lt. u1) goto 10 gamma2 = x*beta return end ************************************************* * * * FUNCTION expon(beta) * * * ************************************************* function expon(beta) implicit real*8 (a-h,o-z) a = 0.0d0 10 u = rand() u0 = u 20 ustar = rand() if(u .lt. ustar) goto 30 u = rand() if(u .lt. ustar) goto 20 a = a+1.0d0 goto 10 30 expon = (a+u0)*beta return end ********************************************************************** * * * FUNCTION gama * * * ********************************************************************** real*8 function gama(x) implicit real*8 (a-h, o-z) b = x a = 1.0d0/(x*(x+1.0d0)) 1 b = b - 1.0d0 if(b.lt.0.0d0) goto 2 a = a*(b+2.0d0) goto 1 2 n=x d=n b=x-d p=(((-0.453010476515162d2*b-0.261028495061646d3)*b- & 0.181779919115643d4)*b-0.495476897267727d4)*b- & 0.135726501110369d5 q = ((((((0.1d1*b-0.161616062977346d2)*b+ & 0.937769173224393d2)*b- & 0.124794725622738d3)*b-0.931178504559348d3)*b+ & 0.344069924134645d4)*b+0.783534880055955d3)*b- & 0.135726501110369d5 gama = a*p/q return end ********************************************************************** * * * SUBROUTINE gc * * * ********************************************************************** subroutine gc (rd,de1,c1) ********************************************************************** * use: compute determinant and inverse of covariance matrix for * * arfima(0,d,0) process i.e. inverse and determinant of * * sigma22(n+nq,n+nq) * * parameters of routine: * * rd *input * scalar. contains current value of fractional * * difference parameter d * * de1 *output* scalar. contains determinant of the matrix sigma22 * * c1 *output* array of length nnq*(nnq+1)/2. contains elements of * * the inverse of sigma22 * * other routines used: gama * ********************************************************************** parameter (nmx=400,npqmx=10) parameter (nmx2=42230) implicit real*8 (a-h,o-z) real*8 c(nmx2) real*8 c1(nmx2) real*8 phi(0:nmx+npqmx) common/nconst/np,nq,ixd,npq,npqd,n,nnq,iter1,nrep1,nloop c de1=0.0d0 if (ixd.eq.0) return c c nnq=n+nq c nd=int((nnq)/2) nd1=int((nnq+1)/2) ********************************************************************* * calculate c(i,j) in order to get inverse matrix * * get the partial linear regression coefficients * * hosking(1981) * * phi(n-1,j)=phi(n-1,j-1)*(n-j)*(j-d-1)/(j)(n-d-j) * ********************************************************************* phi(0)=-1. phi(nnq)=0. phi(1)=(nnq-1)*rd/((nnq-1)*1.-rd) do 2000 j=2,nnq-1 xd1=(j-1)*1.-rd xd2=(nnq-j)*1.-rd phi(j)=phi(j-1)*(nnq-j)*xd1/(j*xd2) 2000 continue c d1=1.-rd d2=1.-2*rd ***************************************** * vt(1)**2=gamma(0)/sigma**2 * ***************************************** vt=gama(d2)/(gama(d1)**2) de1=dlog(vt) ***************************************** * vt(i+1)=vt(i)(1-phi(i)**2) * ***************************************** do 2200 i=1,nnq-1 vt=vt*(1.-(rd/(i*1.-rd))**2) 2200 de1=de1+dlog(vt) c do 2500 i=1,nd1 do 2500 j=i,nnq-i+1 call gctran (i,j,ij,nnq) 2500 c(ij)=(phi(i-1)*phi(j-1)-phi(nnq-i+1)*phi(nnq-j+1))/vt c c do 105 j=1,nnq 105 c1(j)=c(j) c do 110 i=2,nd1 do 110 j=i,nnq-i+1 ij = (nnq-i+2)*(i-1)+j-i+1 110 c1(ij)=c1((nnq-i+3)*(i-2)+j-i+1) + c(ij) c return end ********************************************************************** * * * SUBROUTINE gcova * * * ********************************************************************** subroutine gcova (cova,rd) ********************************************************************** * use: compute the acf of arfima(0,d,0) process * * (hosking, biometrika, 1981) * * parameters of routine: * * cova *output* array of length nw. contains the acf of (0,d,0) * * rd *input * scalar. contains current value of fractional * * difference parameter d * * other routines used: gama * ********************************************************************** parameter (nw=500) implicit real*8 (a-h,o-z) real*8 cova(0:nw) common/nconst/np,nq,ixd,npq,npqd,n,nnq,iter1,nrep1,nloop c d1=1.0d0-rd d2=1.0d0-2.0d0*rd cova(0)=gama(d2)/(gama(d1)**2) do 101 i=1,nw 101 cova(i)=cova(i-1)*(i*1.-1.+rd)/(i*1.-rd) c return end ***************************************** * * * SUBROUTINE gctran * * * ***************************************** subroutine gctran (i,j,ij,n) ********************************************************************** * vectorize a double symmetric matrix * ********************************************************************** if (i.gt.j) then if ((i+j).gt.(n+1)) then ij = (i+1)*(n-i)-j+i+1 else ij = (n-j+2)*(j-1)-j+i+1 endif else if ((i+j).gt.(n+1)) then ij = (j+1)*(n-j)+j-i+1 else ij = (n-i+2)*(i-1)+j-i+1 endif endif c return end ********************************************************************** * * * SUBROUTINE gensig * * * ********************************************************************** subroutine gensig (a1,a2,vari) ********************************************************************** * use: generate sample of white noise standard deviation * * parameters of routine: * * a1 *input * scalar, sum of squares term in the g1 function * * a2 *input * scalar, sum of squares term in the g2 function * * vari *output* scalar, generated value of white noise variance * * other routines used: gammad * ********************************************************************** parameter (nmx=400,npqmx=10) implicit real*8 (a-h,o-z) common/nconst/np,nq,ixd,npq,npqd,n,nnq,iter1,nrep1,nloop c a=(n+npq)/2. c if (np.eq.0) then b=2./a1 else b=2./(a1+a2) endif c rg=gammad(a,b) vari=1./rg return end ********************************************************************** * * * SUBROUTINE gsy * * * ********************************************************************** subroutine gsy (cova,runphi,runthe,s,covy) ********************************************************************** * use: form the acf of arfima(p,d,q) process and elements of eta(k) * * (pai and ravishanker, 1994). these are used in computing the * * sigma11(np,np) and sigma12(np,n+nq) matrices that appear in * * the g1 and g2 functions * * parameters of routine: * * cova *input * array of len(0:nw) contains acf of fi(0,d,0) * * runphi *input * array of len(0:np) contains current phi value * * runthe *input * array of len(0:nq) contains current theta value * * s *output* array of len(1-np-n:nq-1). eta(k) elements * * where s(k)=eta(k)=cov(y(t),a(t-k)) * * covy *output* array of len(0:np-1). acf of arfima(p,d,q) * * other routines used: armaco * ********************************************************************** parameter (nmx=400,npqmx=10) parameter (nw=500,nco=(npqmx+1)**2) implicit real*8 (a-h,o-z) real*8 runphi(0:npqmx),runthe(0:npqmx) real*8 cov(nw),w(nco) real*8 cova(0:nw),covy(0:npqmx-1) real*8 s(1-npqmx-nmx:npqmx-1) real*8 psi(1-npqmx:nw) common/nconst/np,nq,ixd,npq,npqd,n,nnq,iter1,nrep1,nloop c maxpq=max0(np,nq)+1 c nnq=n+nq c ********************************************************************** * form the psi weights(box&jenkins,p134) * ********************************************************************** if (np.ge.2) then do 4410 i=1-np,-1 4410 psi(i)=0. endif c psi(0)=1. c nk=max0(np-1,nq) do 4100 j=1,nk psi(j)=0.0d0 do 4110 i=1,np 4110 psi(j)=psi(j)+runphi(i)*psi(j-i) 4100 if (j.le.nq) psi(j)=psi(j)-runthe(j) ********************************************************************** * short memory arma(p,q) models * ********************************************************************** if (ixd.eq.0) then do 210 i=1-np,nq-1 210 s(i)=psi(i) call armaco (cov,maxpq,np,nq,runphi,runthe,1.0d0,w,nco) c do 220 i=1,np 220 covy(i-1)=cov(i) return endif ********************************************************************** * long memory arfima(p,d,q) models * ********************************************************************** do 200 k=nk+1,nw psi(k)=0.0d0 do 544 l=1,np 544 psi(k)=psi(k)+runphi(l)*psi(k-l) 200 continue c do 310 l=1-np-n,-n,1 s(l)=cova(-l) do 310 i=1,nw+l 310 s(l)=s(l)+cova(iabs(i-l))*psi(i) c do 320 l=1-n,nq-1 s(l)=0. do 330 i=1,np 330 s(l)=s(l)+runphi(i)*s(l-i) do 340 i=0,nq 340 s(l)=s(l)-runthe(i)*cova(iabs(l-i)) 320 continue c call armaco (cov,nw,np,nq,runphi,runthe,1.0d0,w,nco) c do 920 k=0,np-1 covy(k)=cov(1)*cova(k) do 910 j=1,nw-np 910 covy(k)=covy(k)+cov(j+1)*(cova(j+k)+cova(iabs(k-j))) 920 continue c return end ********************************************************************** * * * SUBROUTINE musig * * * ********************************************************************** subroutine musig (phi,the,y0,e0,data,runeps) ********************************************************************** * use:form the error terms: runeps from arfima(0,d,0) process * * of length 1-nq:n * * parameters of routine: * * phi *input * array of length 0:np * * the *input * array of length 0:nq * * y0 *input * array of length 1-np:0 * * e0 *input * array of length 1-nq:0 * * data *input * array of length n * * runeps *output* array of length 1-nq:n * ********************************************************************** parameter (npqmx=10,nmx=400) implicit real*8 (a-h,o-z) real*8 phi(0:npqmx),y0(1-npqmx:0) real*8 data(nmx),runeps(1-npqmx:nmx) real*8 the(0:npqmx),runy(1-npqmx:nmx),e0(1-npqmx:0) common/nconst/np,nq,ixd,npq,npqd,n,nnq,iter1,nrep1,nloop c if (np.ne.0) then do 71 i=1-np,0 71 runy(i)=y0(i) endif c if (nq.ne.0) then do 73 i=1-nq,0 73 runeps(i)=e0(i) endif c do 72 i=1,n 72 runy(i)=data(i) c do 8 mt=1,n runeps(mt)=runy(mt) if (np.ne.0) then do 9 i=1,np 9 runeps(mt)=runeps(mt)-phi(i)*runy(mt-i) endif if (nq.ne.0) then do 10 i=1,nq 10 runeps(mt)=runeps(mt)+the(i)*runeps(mt-i) endif 8 continue c return end ********************************************************************** * * * SUBROUTINE perturb * * * ********************************************************************** subroutine perturb (mphi1,mthe1,y01,eps01,sigma1, / mean1,rd1,rsig,omgch,ap1,nap1) ***************************************************************** * use: perturb the initial values to obtain * * starting values for nrep1 parallel gibbs chains. * * gelman and rubin idea of over-dispersed distribution is used * * parameters of routine: * * mphi1 *input/output* array of dim. (iter1,nrep1,np) * * mthe1 *input/output* array of dim. (iter1,nrep1,nq) * * y01 *input/output* array of dim. (iter1,nrep1,1-np:0) * * eps01 *input/output* array of dim. (iter1,nrep1,1-nq:0) * * sigma1 *input/output* array of dim. (iter1,nrep1) * * mean1 *input/output* array of dim. (iter1,nrep1) * * rd1 *input/output* array of dim. (iter1,nrep1) * * rsig *input * array of dim. (npqd,npqd). contains * * Cholasky dec. of var(phi,theta,d) * * omgch *input * array of len. (npq*(npq+1)/2) * * Cholasky dec. of var(eps0,yo) * * other routines used: unormd, rand, troot * ***************************************************************** c parameter (nrepx=10,iterx=9000,npqmx=10) implicit real*8 (a-h,o-z) real*8 ap1(nap1) real*8 mphi1(iterx,nrepx,npqmx),mthe1(iterx,nrepx,npqmx) real*8 y01(iterx,nrepx,1-npqmx:0) real*8 eps01(iterx,nrepx,1-npqmx:0) real*8 sigma1(iterx,nrepx),mean1(iterx,nrepx) real*8 rd1(iterx,nrepx) real*8 tempp(0:npqmx),tempq(0:npqmx) real*8 omgch(npqmx*(npqmx+1)/2) real*8 r(npqmx),rsig(npqmx*(npqmx+1)/2) common/nconst/np,nq,ixd,npq,npqd,n,nnq,iter1,nrep1,nloop data overdis/1.0d0/ c if (ixd.eq.1) rd1(1,1)=ap1(npqd) c mean1(1,1)=ap1(npqd+1) sigma1(1,1)=ap1(npqd+3) c if (np.ne.0) then do 111 i=1,np y01(1,1,i-np)=ap1(npqd+1) 111 mphi1(1,1,i)=ap1(i) endif c if (nq.ne.0) then do 112 i=1,nq eps01(1,1,i-nq)=0.0d0 112 mthe1(1,1,i)=ap1(i+np) endif c if (nrep1.eq.1) return c tempp(0)=-1.0d0 tempq(0)=-1.0d0 c do 10 i=2,nrep1 c if (npq.ne.0) then call dmvnor (omgch,r,npq) if (nq.ne.0) then do 21 j=1,nq eps01(1,i,j-nq)=eps01(1,1,j-nq)+r(j+np)*overdis 21 continue endif c if (np.ne.0) then do 20 j=1,np y01(1,i,j-np)=y01(1,1,j-np)+r(j)*overdis 20 continue endif endif c 2000 call dmvnor (rsig,r,npqd) c if (np.ne.0) then do 22 j=1,np tempp(j)=mphi1(1,1,j)+r(j)*overdis 22 continue call troot (tempp,np,indicat) if (indicat.eq.0) goto 2000 endif c if (nq.ne.0) then do 23 j=1,nq tempq(j)=mthe1(1,1,j)+r(j+np)*overdis 23 continue call troot (tempq,nq,indicat) if (indicat.eq.0) goto 2000 endif c if (ixd.eq.1) then rd1(1,i)=rd1(1,1)+r(npqd)*overdis if (dabs(rd1(1,i)).ge.0.5d0) goto 2000 c call dblepr ("rd1",3,rd1(1,i),1) endif c if (np.ne.0) then do 1210 j=1,np mphi1(1,i,j)=tempp(j) 1210 continue endif c if (nq.ne.0) then do 310 j=1,nq mthe1(1,i,j)=tempq(j) 310 continue endif c mean1(1,i)=mean1(1,1)+unormd()*overdis*ap1(npqd+2) sigma1(1,i)=sigma1(1,1)*(1.0d0+rand()) c 10 continue c return end function rand() ************************************************* * * * This function returns a random uniform * * number (0 - 1) based multiplicative * * congruential method. * * * * Note: The initial seed is 123456789 * * (it must be odd value). * * * * The period is 2^31(=2147483648). * * * ************************************************* implicit real*8 (a-h, o-z) real*8 invm parameter( lamda = 48828125, mask = 2**30+(2**30-1) ) parameter( invm = 0.5d0**31 ) data is /0/ save iz, is c if(is.eq.0) then is = 1 ix = 123456789 iz = ix end if c iz = jiand(lamda * iz, mask) rand = iz*invm return end ********************************************************************** * revise this sub !!!!!! * * SUBROUTINE troot * * * ********************************************************************** subroutine troot (runbet,nbet,indicat) ********************************************************************** * test whether a model is stationary for given coefficients * * nbet=np+1 or nq+1 * ********************************************************************** parameter (npqmx=10) implicit real*8 (a-h,o-z) real*8 runbet(0:nbet),a(npqmx),b(npqmx) integer indicat c do 110 i=0,nbet 110 a(i+1)=runbet(i) a(1)=-1. do 2 k=nbet,1,-1 do 3 j=1,k b(j)=a(1)*a(j)-a(k+1)*a(k+2-j) 3 continue if (b(1).le.0.0)go to 200 do 4 j=1,k a(j)=b(j) 4 continue 2 continue indicat=1 go to 220 200 indicat=0 220 return end ********************************************************************** * * * SUBROUTINE updatemu * * * ********************************************************************** subroutine updatemu (fu,c1,de1,runmn,runphi,runthe,stdmu / ,data,runy0,runeps0,runeps,vari,rd,cova,covy,s) ********************************************************************** * use: generate sample of mu by metropolis algorithm * * parameters of routine: * * fu *input/output* scalar. log of posterior density * * c1 *input * array of dimension (n+nq,n+nq). elements * * of inverse of sigma22 * * de1 *input * scalar. determinant in the g1 function * * runmn *input/output* scalar. generated sample of mean mu * * runphi *input * array of length 0:np. * * runthe *input * array of length 0:nq * * stdmu *input * scalar. estimated standard deviation for the * * proposal density for sampling mu * * data *input/output* array of length n : mean subtracted data * * runy0 *input/output* array of length 1-np:0 * * runeps0 *input * array of length 1-nq:0 * * runeps *input * array of length 1-nq:n * * vari *input * scalar. white moise variance * * rd *input * scalar. fractional difference parameter d * c cova *input * array of length 0:nw. acf from fi(d) process * c covy *input * array of length 0:np-1. acf from arfima(p,d,q) * c process * c s *input * array of length 1-np-n:nq-1. contains eta(k) * c other routines used: unormd,rand,musig,g1,g2,fcal * ********************************************************************** parameter (npqmx=10) parameter (nmx=400,nw=500,nloopx=20) parameter (nmx2=42230) implicit real*8 (a-h,o-z) real*8 runphi(0:npqmx) real*8 runthe(0:npqmx) real*8 rn(nloopx) real*8 data(nmx),runeps(1-npqmx:nmx) real*8 runepst(1-npqmx:nmx) real*8 runy0(1-npqmx:0),runeps0(1-npqmx:0) real*8 c1(nmx2) real*8 cova(0:nw),covy(0:npqmx-1) real*8 s(1-npqmx-nmx:npqmx-1) common/nconst/np,nq,ixd,npq,npqd,n,nnq,iter1,nrep1,nloop c do 1000 loop=1,nloop rn(loop)=unormd()*stdmu runmn=runmn+rn(loop) c do 5 i=1,n data(i)=data(i)-rn(loop) 5 continue c if (np.ne.0) then do 15 i=1,np 15 runy0(i-np)=runy0(i-np)-rn(loop) endif call musig (runphi,runthe,runy0,runeps0,data,runepst) call g1 (runepst,c1,a1) if (np.ne.0) call g2 (c1,runepst,runy0,covy,s,rd,a2,de2) call fcal (fv,a1,a2,de1,de2,vari) c write(6,*)'updatemu a1 a2 de1 de2',a1,a2,de1,de2 c fuv=fv-fu ulg=alog(rand()) c write(6,*)'updatemu fu,fv,fuv,ulg',fu,fv,fuv,ulg if (ulg.lt.fuv) then fu=fv do 100 i=1-nq,n 100 runeps(i)=runepst(i) else runmn=runmn-rn(loop) do 8 i=1,n 8 data(i)=data(i)+rn(loop) if (np.ne.0) then do 18 i=1,np 18 runy0(i-np)=runy0(i-np)+rn(loop) endif endif c 1000 continue return end ********************************************************************** * * * SUBROUTINE update * * * ********************************************************************** subroutine update (fu,c1,de1,rd,runphi,runthe,runy0,runeps0 / ,vari,a1,a2,rsig,runeps,data,cova,covy,s) ********************************************************************** * use: generate sample for (PHI,THETA,D) vectors * * usig metropolis algorithm * * parameters of routine: * * fu *input/output* scalar. log of posterior density * * c1 *input/output* array of dimension (n+nq,n+nq). contains * * elements of inverse of sigma22 * * de1 *input/output * scalar. contains determinant term in g1 * * function * * rd *input/output* scalar. generated value of fractional * * difference parameter d * * runphi *input/output* array of length 0:np. contains generated * * sample for phi vector * * runthe *input/output* array of length 0:nq. * * runy0 *input * array of length 1-np:0 * * runeps0 *input * array of length 1-nq:0 * * vari *input * scalar. whtie noise variance * * a1 *input/output* scalar. sum of squares term in g1 function * * a2 *input/output* scalar. sum of squares term in g2 function * * rsig *input* array of dimension (npqd*(npqd+1)/2) * * where npqd=npq+ixd * * contains the Cholesky decoposition of * * VAR(PHI,THETA,d) * * runeps *input/output* array of length 1-nq:n * * data *input * array of length n * * cova *input * array of length 0:nw. acf from * * arfima(0,d,0) process * * covy *input * array of length 0:np-1. acf from * * arfima(p,d,q) process * * s *input * array of length 1-np-n:nq-1. eta(k) * * other routines used: rand,dmvnor,troot,musig,g1,gsy * * g2,fcal * ********************************************************************** parameter (npqmx=10,nmx=400,nloopx=20) parameter (nmx2=42230) parameter (nw=500) implicit real*8 (a-h,o-z) real*8 runphi(0:npqmx) real*8 runthe(0:npqmx) real*8 tempp(npqmx) real*8 tempq(npqmx) real*8 runepst(1-npqmx:nmx) real*8 data(nmx),runeps(1-npqmx:nmx) real*8 runy0(1-npqmx:0),runeps0(1-npqmx:0) real*8 rsig(npqmx*(npqmx+1)/2),r(npqmx),vari real*8 c1(nmx2) real*8 c1t(nmx2) real*8 cova(0:nw),covy(0:npqmx-1) real*8 s(1-npqmx-nmx:npqmx-1) real*8 covat(0:nw) real*8 covyt(0:npqmx-1) real*8 st(1-npqmx-nmx:npqmx-1) common/nconst/np,nq,ixd,npq,npqd,n,nnq,iter1,nrep1,nloop c ********************************************************************** inmx2=int(((nnq+1)**2)/4) loop=0 c 2000 if (np.ne.0) then do 1 i=1,np 1 tempp(i)=runphi(i) endif c if (nq.ne.0) then do 2 i=1,nq 2 tempq(i)=runthe(i) endif c if (ixd.eq.1) tempd=rd c 1000 call dmvnor (rsig,r,npqd) c if (np.ne.0) then do 5 i=1,np runphi(i)=runphi(i)+r(i) 5 continue ********************************* * check for stationarity * ********************************* call troot (runphi,np,indicat) if (indicat.eq.0) then do 105 i=1,np runphi(i)=runphi(i)-r(i) 105 continue goto 1000 endif endif c if (nq.ne.0) then do 7 i=1,nq runthe(i)=runthe(i)+r(i+np) 7 continue ********************************* * check for invertibility * ********************************* call troot (runthe,nq,indicat) if (indicat.eq.0) then do 106 i=1,nq runthe(i)=runthe(i)-r(i+np) 106 continue goto 1000 endif endif c if (ixd.eq.1) then rd=rd+r(npqd) c ******************************************************* * check for stationarity and invertibility * ******************************************************* if (dabs(rd).ge.0.5d0) then rd=rd-r(npqd) goto 1000 endif endif ********************************************************* * satisfies stationarity and invertibility * ********************************************************* loop=loop+1 call musig (runphi,runthe,runy0,runeps0,data,runepst) call gc (rd,de1t,c1t) call g1 (runepst,c1t,a1t) if (np.ne.0) then if (ixd.eq.1) call gcova (covat,rd) call gsy (covat,runphi,runthe,st,covyt) call g2 (c1t,runepst,runy0,covyt,st,rd,a2t,de2t) endif call fcal (fv,a1t,a2t,de1t,de2t,vari) c write(6,*)'update a1 a2 de1 de2',a1t,a2t,de1t,de2t c fuv=fv-fu c ulg=alog(rand()) c write(6,*)'update fu,fv,fuv,ulg',fu,fv,fuv,ulg if (ulg.lt.fuv) then fu=fv do 200 i=1-nq,n 200 runeps(i)=runepst(i) a1=a1t de1=de1t do 100 ij=1,inmx2 100 c1(ij)=c1t(ij) if (np.ne.0) then a2=a2t de2=de2t if (ixd.eq.1) then do 250 i=0,nw 250 cova(i)=covat(i) endif do 300 i=0,np-1 300 covy(i)=covyt(i) do 400 i=1-np-n,nq-1 400 s(i)=st(i) endif else if (np.ne.0) then do 205 i=1,np runphi(i)=tempp(i) 205 continue endif if (nq.ne.0) then do 207 i=1,nq runthe(i)=tempq(i) 207 continue endif if (ixd.eq.1) rd=tempd end if c if (loop.lt.nloop) goto 2000 c return end ********************************************************************** * * * SUBROUTINE upeps * * * ********************************************************************** subroutine upeps (fu,c1,de1,runphi,runthe,runy0,runeps0 / ,omgch,data,runeps,vari,rd,cova,covy,s) ********************************************************************** * use: generate sample for latent variables y0 and eps0 using * * metropolis algorithm * * parameters of routine: * * fu *input/output* scalar. log of posterior density * * c1 *input* array of dimension (n+nq,n+nq). contains * * elements of inverse of sigma22 * * de1 *input * scalar. contains determinant term in g1 function* * runphi *input * array of length 0:np. * * runthe *input * array of length 0:nq. * * runy0 *input/output* array of length 1-np:0. * * generated sample for * * y0 vector (latent variables) * * runeps0 *input/output* array of length 1-nq:0. * * generated sample for * * eps0 vector(latent variables) * * omgch *input * array of dimension (npq*(npq+1)/2). estimated * * Cholesky decoposition of * * for sampling (y0,eps0) * * data *input * array of length n * * runeps *input/output* array of length 1-nq:n * * vari *input * scalar. whtie noise variance * * rd *input * scalar. fractional difference parameter d * * cova *input * array of length 0:nw. acf from * * arfima(0,d,0) process * * covy *input * array of length 0:np-1. acf from * * arfima(p,d,q) process * * s *input * array of length 1-np-n:nq-1. eta(k) * * other routines used: rand,dmvnor,musig,g1,g2,fcal * ********************************************************************** parameter (npqmx=10,nmx=400,nloopx=20) parameter (nmx2=42230) parameter (nw=500) implicit real*8 (a-h,o-z) real*8 r(npqmx) real*8 runphi(0:npqmx),runthe(0:npqmx) real*8 omgch(npqmx*(npqmx+1)/2) real*8 data(nmx),runeps(1-npqmx:nmx) real*8 runepst(1-npqmx:nmx) real*8 runy0(1-npqmx:0),runeps0(1-npqmx:0) real*8 c1(nmx2) real*8 cova(0:nw),covy(0:npqmx-1) real*8 s(1-npqmx-nmx:npqmx-1) common/nconst/np,nq,ixd,npq,npqd,n,nnq,iter1,nrep1,nloop c do 50 loop=1,nloop call dmvnor (omgch,r,npq) if (nq.ne.0) then do 52 i=1,nq 52 runeps0(i-nq)=runeps0(i-nq)+r(i+np) endif c if (np.ne.0) then do 54 i=1,np 54 runy0(i-np)=runy0(i-np)+r(i) endif c call musig (runphi,runthe,runy0,runeps0,data,runepst) call g1 (runepst,c1,a1) if (np.ne.0) call g2 (c1,runepst,runy0,covy,s,rd,a2,de2) call fcal (fv,a1,a2,de1,de2,vari) c write(6,*)'upeps a1 a2 de1 de2',a1,a2,de1,de2 fuv=fv-fu ulg=alog(rand()) c write(6,*)'upeps fu,fv,fuv,ulg',fu,fv,fuv,ulg c if (ulg.lt.fuv) then fu=fv do 100 i=1-nq,n 100 runeps(i)=runepst(i) else if (nq.ne.0) then do 152 i=1,nq 152 runeps0(i-nq)=runeps0(i-nq)-r(i+np) endif if (np.ne.0) then do 154 i=1,np 154 runy0(i-np)=runy0(i-np)-r(i) endif endif c 50 continue c return end ********************************************************************** * * * SUBROUTINE xomega * * * **************************************************************************** subroutine xomega (omgch,ap1,nap1) ********************************************************************** * use: form the covariance matrix of the latent variables eps0,y0 * * by a method similar to that used in newbold(1974) for arma * * models and then convert to its Cholesky decoposition * * parameters of routine: * * omega *output* array of dim. (npq,npq). Cholesky decoposition * * of VAR(eps0,y0) * * phi *input * array of length 0:np * * the *input * array of length 0:nq * * rd *input * scalar fractional difference parameter d * * other routines used: gama, armaco * ********************************************************************** parameter (npqmx=10,nmx=400,nw=500) parameter (nco=(npqmx+1)**2) implicit real*8 (a-h,o-z) real*8 psi(1-npqmx:nw),omgch(npqmx*(npqmx+1)/2) real*8 the(0:npqmx),phi(0:npqmx) real*8 cov(nw),w(nco),omega(npqmx,npqmx) real*8 cova(0:nw),covy(0:npqmx-1) real*8 s(1-npqmx-nmx:npqmx-1),ap1(nap1) common/nconst/np,nq,ixd,npq,npqd,n,nnq,iter1,nrep1,nloop c maxpq=max0(np,nq)+1 c if (np.ne.0) then do 20 i=1,np 20 phi(i)=ap1(i) endif c if (nq.ne.0) then do 30 i=1,nq 30 the(i)=ap1(np+i) endif c if (ixd.eq.1) then rd=ap1(np+nq+ixd) d1=1.0d0-rd d2=1.0d0-2.0d0*rd cova(0)=gama(d2)/(gama(d1)**2) endif ********************************************************************** * arfima(p,0,0) models * ********************************************************************** if (nq.eq.0 .and. ixd.eq.0) then call armaco (cov,maxpq,np,0,phi,the,1.0d0,w,nco) do 630 i=1,np do 630 j=1,np 630 omega(i,j)=cov(iabs(i-j)+1) endif ********************************************************************** * arfima(p,d,0) models * ********************************************************************** if (nq.eq.0 .and. ixd.eq.1) then call armaco(cov,nw,np,0,phi,the,1.0d0,w,nco) do 100 i=1,nw 100 cova(i)=cova(i-1)*(float(i)-1.0d0+rd)/(float(i)-rd) c do 120 k=0,np-1 covy(k)=cov(1)*cova(k) do 1001 j=1,nw-np 1001 covy(k)=covy(k)+cov(j+1)*(cova(j+k)+cova(iabs(k-j))) 120 continue c do 530 i=1,np do 530 j=1,np omega(i,j)=covy(iabs(i-j)) 530 continue endif ********************************************************************** * arfima(0,0,q) models * ********************************************************************** if (np.eq.0 .and. ixd.eq.0) then do 610 i=1,nq do 610 j=1,nq if (i.eq.j) then omega(i,j)=1.0d0 else omega(i,j)=0.0d0 endif 610 continue endif ********************************************************************** * arfima(0,d,q) models * ********************************************************************** if (np.eq.0 .and. ixd.eq.1) then do 130 i=1,nq 130 cova(i)=cova(i-1)*(float(i)-1.0d0+rd)/(float(i)-rd) c do 510 i=1,nq do 510 j=1,nq 510 omega(i,j)=cova(iabs(i-j)) endif ********************************************************************** * psi weights for arfima(p,0,q) or arfina(p,d,q) models * * see box & jenkins,p134 * ********************************************************************** if (np.ne.0 .and. nq.ne.0) then if (np.ge.2) then do 4410 i=1-np,-1 4410 psi(i)=0.0d0 endif psi(0)=1.0d0 c nk=max0(np-1,nq) do 4100 j=1,nk psi(j)=0. do 4110 i=1,np 4110 psi(j)=psi(j)+phi(i)*psi(j-i) 4100 if (j.le.nq) psi(j)=psi(j)-the(j) endif ********************************************************************** * arfima(p,0,q) models * ********************************************************************** if (np.ne.0 .and. ixd.eq.0 .and. nq.ne.0) then call armaco(cov,maxpq,np,nq,phi,the,1.0d0,w,nco) do 635 i=1,np do 635 j=1,np omega(i,j)=cov(iabs(i-j)+1) 635 continue c do 640 i=1,nq do 640 j=1,np 640 omega(i,np+j)=psi(j-np-i+nq) c640 omega(np+j,i)=omega(i,np+j) c do 625 i=1,nq do 625 j=1,nq if (i.eq.j) then omega(np+i,np+j)=1.0d0 else omega(np+i,np+j)=0.0d0 endif 625 continue endif ********************************************************************** * arfima(p,d,q) models * ********************************************************************** if (np.ne.0 .and. ixd.eq.1 .and. nq.ne.0) then do 4500 k=nk+1,nw psi(k)=0. do 4510 l=1,np 4510 psi(k)=psi(k)+phi(l)*psi(k-l) 4500 continue c call armaco(cov,nw,np,nq,phi,the,1.0d0,w,nco) do 140 i=1,nw 140 cova(i)=cova(i-1)*(float(i)-1.0d0+rd)/(float(i)-rd) c do 150 k=0,np-1 covy(k)=cov(1)*cova(k) do 2001 j=1,nw-np 2001 covy(k)=covy(k)+cov(j+1)*(cova(j+k)+cova(iabs(k-j))) 150 continue c do 525 i=1,nq do 525 j=1,nq 525 omega(np+i,np+j)=cova(iabs(i-j)) c do 542 l=1-np,nq-1 s(l)=cova(iabs(l)) do 542 k=1,nw-np+1 542 s(l)=s(l)+psi(k)*cova(iabs(k-l)) c do 540 i=1,nq do 540 j=1,np 540 omega(i,np+j)=s(j-np-i+nq) c540 omega(np+j,i)=omega(i,np+j) c do 535 i=1,np do 535 j=1,np omega(i,j)=covy(iabs(i-j)) 535 continue endif c ij=0 do 700 j=1,np+nq do 700 i=1,j ij=ij+1 700 omgch(ij)=omega(i,j)*ap1(nap1)**2 ***************************************************************** call chod (omgch,npq) ***************************************************************** c return end function unormd() ***************************************** * * * This function generates a random * * normal number with mean av=0.0 and * * standard deviation sd=1.0. * * * * * * av ...... mean * * sd ...... standart deviation * * * * * * Reference RIPLEY, B.D. (1987) * * "STOCHASTIC SIMULATION" * * * * * ***************************************** implicit real*8 (a-h, o-z) save ir, an data ir/0/ if(ir.eq.0) then 10 z1 = rand() z2 = rand() v1 = 2.0d0*z1 - 1.0d0 v2 = 2.0d0*z2 - 1.0d0 w = v1*v1 + v2*v2 if(w.gt.1.0d0) goto 10 e = dsqrt( (-2.0d0*dlog(w))/w ) an = v1*e ir = 1 unormd = v2*e else ir = 0 unormd = an end if return end ********************************************************************** * FUNCTION XSMINV ********************************************************************** real*8 function xsminv (y9,n) parameter (npqmx=10) implicit real*8 (a-h,o-z) real*8 x(210),y9(npqmx,npqmx) ********************************************************************** c* * c* purposes: function xsminv inverts a given (n*n)-dimensioned * c* symmetric positive definite matrix, using the * c* cholesky's upper triangular factorization method. * c* * c* usage : det=xsminv(y,n) * c* * c* * c* parameters on entry on return * c* * c* y : symmetric matrix inversed symmtric matrix * c* * c* n : number of rows unchanged * c* * c* remark : error glag xsminv is set to the following values: * c* * c* xsminv > 0.0 : no error - matrix's determination * c* xsminv = 0.0 : singular matrix * c* xsminv = -1.0 : wrong parameter n on entry * c* xsminv = -5.0 : matrix not definite positive * c* * ********************************************************************** c do 1 j=1,n do 1 i=1,j 1 x((j-1)*j/2+i)=y9(i,j) c det=1.0 if (n.le.0) then xsminv=-1.0 return endif c c compute upper-triangular factorization matrix t, c such as x=t'*t, by the cholesky's square root method. c do 1000 i=1,n max=(i*i+i)/2 ijk=max do 1000 j=i,n y=0.0 if (i.gt.1) then do 100 k=1,i-1 l=max-k m=ijk-k 100 y=y+x(l)*x(m) endif y=x(ijk)-y if (j.eq.i) then if (y.lt.0.0) then xsminv=-5.0 return elseif (y.eq.0.0) then xsminv=0.0 return endif z=dsqrt(y) x(max)=z det=det*x(max) z=1.0/z else x(ijk)=y*z endif 1000 ijk=ijk+j c c invert factorized upper-triangular matrix t. c ijk=max do 2000 i=1,n y=1.0/x(max) x(max)=y min=n if (i.gt.1) then klm=ijk do 200 j=1,i-1 z=0.0 min=min-1 l=max m=klm do 20 k=n-i+1,min l=l+k m=m+1 20 z=z+x(l)*x(m) x(klm)=-y*z 200 klm=klm-min endif max=max-min 2000 ijk=ijk-1 c c compute (x)(inv) by means of (t)(inv): c (x)(inv)=(t)(inv)*((t)(inv))' c do 3000 i=1,n max=(i*i+i)/2 ijk=max do 3000 j=i,n z=0.0 m=ijk do 300 k=j,n l=m+j-i z=z+x(l)*x(m) 300 m=m+k x(ijk)=z 3000 ijk=ijk+j xsminv=det*det c c do 9 j=1,n do 9 i=1,j y9(i,j)=x((j-1)*j/2+i) 9 y9(j,i)=y9(i,j) c if (xsminv.eq.0.0) then call dblepr (" 6 singular matrix, xsmins=",50,xsminv,1) endif ********12345678901234567890123456789012345678901234567890" c return end c ********************************************************************** * SUBROUTINE armaco * ********************************************************************** subroutine armaco(cov,ncov,ip,iq,phi,theta,wnvar,w,nw) ********************************************************************** c covariance or correlation function of an arma process c c parameters of routine: c cov *output array of length at least ncov. if wnvar.gt.0, then o c exit c(1) contains the variance of the arma process c and c(2)...c(ncov) contain the covariances at lags c l...(ncov-1). if wnvar.le.0, then on exit c(1)=1 and c c(2)...c(ncov) contain the correlations at lags c 1...(ncov-1). c input * number of covariances to be found. must be greater c ncov * than max(ip,iq). c ip * input* number of autoregressive parameters c iq * input* number of moving-average parameters c phi * input* array of length ip. contains the autoregressive c parameters. c theta * input* array of length iq. contains the moving-average c parameters. c wnvar * input* if positive, the white-noise variance of the arma c process. if negative or zero, correlations rather c than covariances wttit. e found. c c c w * local* work array of lengtii nw c nw * input* size of work array. must be at least: c (if ip.eq.o) 1; c (if ip.eq.1) max(1,iq*2); c (if ip.ge.2) (m+1)**2, where m=max(ip,iq). c c method used is that of mcleod (appl.statist. 1975,1977) except that c the signs of the a and b coefficients are reversed. c no check is made that the arma process is stationary or invertible. c other routines used: determ,lineq ********************************************************************** implicit real*8 (a-h,o-z) parameter (npqmx=10,nwx=500) real*8 cov(nwx),phi(0:npqmx),theta(0:npqmx),w(nw) do 10 k=1,ncov 10 cov(k)=0. if(ncov.le.ip.or.ncov.le.iq)goto 1000 if(ip-1)20,15,110 15 if(iq-1)70,80,90 20 if(iq.gt.0)goto 30 c c arma(0,0) c cov(1)=wnvar if(wnvar.le.0.)cov(1)=1. return c c arma(0,q) (pure moving average) c 30 sum=1. do 35 i=1,iq 35 sum=sum+theta(i)**2 c0=wnvar if(wnvar.le.0.)c0=1./sum cov(1)=sum*c0 if(iq.eq.1)goto 60 kk=1 do 50 k=2,iq sum=-theta(kk) do 40 i=k,iq 40 sum=sum+theta(i)*theta(i-kk) cov(k)=sum*c0 50 kk=k 60 cov(iq+1)=-theta(iq)*c0 return c c arma(1,o) c 70 p=1-phi(1)**2 cov(1)=wnvar/p if(wnvar.le.0.)cov(1)=1. do 75 k=2,ncov 75 cov(k)=cov(k-1)*phi(1) return c c arma(1,1) c 80 p=1-phi(1)**2 cov(1)=(1.+theta(1)**2-2.*phi(1)*theta(1))/p c0=wnvar if(wnvar.le.0.)c0=1./cov(1) cov (1)=cov(1)*c0 cov(2)=(1-phi(1)*theta(1))*(phi(1)-theta(1))*c0/p if(ncov.eq.2)return do 85 k=3,ncov 85 cov(k)=cov(k-1)*phi(1) return c c arma(1,q) c 90 if(nw.lt.iq+iq)goto 1010 iqq=iq+1 p=phi(1) c c find psi weights c w(1)=1. jj=1 do 92 j=2,iqq w(j)=p*w(jj)-theta(jj) 92 jj=j do 94 j=iq+2,iq+iq w(j)=p*w(jj) 94 jj=j c c find covariances c pp=w(iqq)**2/(1-p*p) ii=0 do 98 i=1,iqq sum=w(i)+pp do 96 j=2,iq 96 sum=sum+w(j)*w(j+ii) cov(i)=sum pp=pp*p 98 ii=i c0=wnvar if(wnvar.le.0.)c0=1./cov(1) do 100 i=1,iqq 100 cov(i)=cov(i)*c0 if(ncov.eq.iqq)return do 102 i=iq+2,ncov cov(i)=p*cov(ii) 102 ii=i c c general case: arma(p,q) (p.ge.2) c 110 irr=max0(ip,iq)+1 if(nw.lt.irr*irr)goto 1010 c c find matrix a c do 120 iw=1,nw 120 w(iw)=0 do 130 j=1,irr iw=j+(j-1)*irr inc=-irr w(iw)=1. do 130 k=1,ip iw=iw+inc if(iw.gt.0)goto 130 iw=j+irr inc=irr 130 w(iw)=w(iw)-phi(k) cov(1)=1. c c find c coefficients c icorig=irr*irr if(iq.eq.0)goto 190 w(icorig+1)=phi(1)-theta(1) if(iq.eq.1)goto 160 do 150 k=2,iq c=-theta(k) do 140 i=1,min0(ip,k) cc=1. if(i.lt.k)cc=w(icorig+k-i) 140 c=c+cc*phi(i) 150 w(icorig+k)=c 160 continue c find b coefficients (stored in array cov) do 180 k=1,iq b=1. if(k.gt.1)b=-theta(k-1) do 170 i=k,iq 170 b=b-theta(i)*w(icorig+i+1-k) 180 cov(k)=b cov(iq+1)=-theta(iq) c c solve a.cov = b c 190 d1=-1. call lineq(w,irr,irr,cov,0,ifail) if(ifail.eq.2)goto 1020 c if(ifail.eq.3) then call intpr (" 6 error in lineq: result may be', inaccurate, ifail=",50,ifail,1) ********12345678901234567890123456789012345678901234567890" endif c c0=wnvar if(wnvar.le.0.)c0=1./cov(1) do 195 i=1,irr 195 cov(i)=cov(i)*c0 if(ncov.eq.irr)return c c find higher-lag covariances from difference equation c do 210 k=irr+1,ncov sum=0. do 200 i=1,ip 200 sum=sum+phi(i)*cov(k-i) 210 cov(k)=sum return c 1000 call intpr (" 6 error in armaco: ncov <= np or ncov <= nq, ncov=",50,ncov,1) return 1010 call intpr (" 6 error in armaco: insufficient workspace ",50,nw,1) return 1020 call intpr (" 6 error in armaco: process may be nonstationary ",50,nw,1) ********12345678901234567890123456789012345678901234567890" cov(1)=1. do 1030 i=2,ncov 1030 cov(i)=0. return end ********************************************************************** * SUBROUTINE dterm * ********************************************************************** subroutine dterm(a,na,n,det,idet,ifail) ********************************************************************** c determinant of a square matrix c c (originally 'subroutine determ': renamed to avoid conflict with c r.f.voss's vp package. used by invmat, lineq.) c c parameters of routine: c a *in/out* 2-dimensional array containing the matrix. on c successful exit, contains the inverse matrix. c na * input* the first dimension of array a, as defined in the c calling program c n * input* the order of the matrix c det *output* mantissa of determinant ) determinant is c idet *output* exponent of determinant ) det*256**idet c ifail *output* error flag. on exit, it is set to: c 0 on successful exit c 1 if na is less than n c 2 if matrix is singular c 3 if matrix is nearly singular c if ifail=2 or 3, the value of the determinant may b c inaccurate. c c c method used is iteration of the formula for the determinant of a c partitioned matrix: c c a b -1 c det ( ) = det(d). det(a - bd c) c c d c c with d scalar. to avoid problems with overflow when dealing with c matrices whose elements are large or small, the determinant is stor c in the form det*e**idet where the base e is 256. the routine is n c particularly accurate and should be converted to extended precision c if greater accuracy is required. ********************************************************************** implicit real*8 (a-h,o-z) real*8 a(na,n) data zero/0d0/ c c e is the base for exponentiation, i.e. determinant is stored c det * e**idet c ei is 1/e c eps controls the test for near-singularity c data e,ei/256d0,0.00390625d0/,eps/1d-12/ if(na.lt.n)goto 1000 amm=a(n,n) if(amm.eq.0d0)goto 1010 det=amm idet=0 10 if(dabs(det).le.e)goto 20 det=det*ei idet=idet+1 goto 10 20 if(dabs(det).ge.ei)goto 30 det=det*e idet=idet-1 goto 20 30 if(n.eq.1)return c mz=n do 80 l=2,n m=mz mz=m-1 do 40 j=1,mz do 40 i=1,mz 40 a(i,j)=a(i,j)-a(i,m)*a(m,j)/amm amm=a(mz,mz) if(amm.eq.zero)goto 1010 det=det*amm 60 if(dabs(det).le.e)goto 70 det=det*ei idet=idet+1 goto 60 70 if(dabs(det).ge.ei)goto 80 det=det*e idet=idet-1 goto 70 80 continue c c test for near-singularity c amm=dabs(a(1,1)) zmin=amm zmax=amm do 90 m=2,n amm=dabs(a(m,m)) if(amm.gt.zmax)zmax=amm if(amm.lt.zmin)zmin=amm 90 continue ifail=0 if(zmin/zmax.ge.eps)return ifail=4 call intpr (" 6 error in cterm: matrix is almost singular, ifail=",50,ifail,1) return c 1000 ifail=1 det=zero idet=0 call intpr (" 6 error in cterm: matrix dimen. incompatible, ifail=",50,ifail,1) return c 1010 ifail=2 det=zero idet=0 call intpr (" 6 error in cterm: matrix is singular(warning) ifail=",50,ifail,1) ********12345678901234567890123456789012345678901234567890" return end ********************************************************************** * SUBROUTINE lineq * ********************************************************************** subroutine lineq (a,na,n,u,iflag,ifail) ********************************************************************** c solves the linear equations a . x = u where a is an n x n matri c and x and u are n-vectors c c parameters of routine: c a *in/out* 2-dimensional array containing the matrix. both c dimensions of a must be .ge. n. the elements of c the array are altered by the routine. c na * input* the first dimension of array a, as defined in the c calling program c n * input* the order of matrix a c u *in/out* array of length n. on entry should contain the vector c u. on successful exit, contains the solution vecto c x. c iflag * input* should usually be set to zero - see below c ifail *output* error flag. on exit, it is set to c 0 on successful exit c 1 if na is less than n c 2 if matrix a is singular c 3 if matrix a is nearly singular c if ifail=2 or 3, the solution may be inaccurate. c c iflag parameter - circumstances in which it may be nonzero. c (1) the determinant of matrix a has previously been calculated usin c routine determ and the elements of array a have not subsequentl c been altered. in this case lineq may be called with iflag=1, c substantial computation being thereby saved. c (2) several sets of equations are to be solved with the same a c matrix but different u vectors. in this case lineq should be c called the first time with iflag=0, and the second and c subsequent times with iflag=1. array a and the variable ifail c should be left unchanged between successive calls to lineq. c again, substantial computation is saved. c c the method used is iteration of the formula for solution of a c partitioned set of linear equations: c c -1 -1 -1 -1 c if a b . x = u then x = (a - b d c) (u - b d v), y = d (v-cx c c d y v c c with d scalar. c c other routines used: dterm ********************************************************************** implicit real*8 (a-h,o-z) real*8 a(na,n),u(n) if(na.lt.n)goto 1000 if(iflag.eq.0)call dterm(a,na,n,det,idet,ifail) if(ifail.eq.2)goto 1010 c if(ifail.eq.3) then call intpr (" 6 error in lineq: result may be', inaccurate, ifail=",50,ifail,1) endif c if(n.eq.1)goto 30 c c -1 c form u - b d v for successively smaller orders c m=n do 20 mm=2,n mz=m-1 r=u(m)/a(m,m) do 10 i=1,mz 10 u(i)=u(i)-r*a(i,m) 20 m=mz c c form solution vector c 30 u(1)=u(1)/a(1,1) if(n.eq.1)return do 50 m=2,n sum=u(m) do 40 i=1,m-1 40 sum=sum-u(i)*a(m,i) 50 u(m)=sum/a(m,m) return c 1000 ifail=1 call intpr (" 6 error in lineq: first dim. of arrau too small, i=",50,ifail,1) return 1010 call intpr (" 6 error in lineq: abandoned, ifail=",50,ifail,1) return end subroutine prx0 (mphi1,mthe1,y01,eps01,mean1,sigma1,rd1) c use:print the samples from gibbs sampler c prints generated samples for each model parameter together with c posterior means and standard deviations c parameters of routine: c mphi1 *input* array of dimension (iter1,nrep1,np) c mthe1 *input* array of dimension (iter1,nrep1,nq) c y01 *input* array of dimension (iter1,nrep1,1-np:0) c eps01 *input* array of dimension (iter1,nrep1,1-nq:0) c mean1 *input* array of dimension (iter1,nrep1) c sigma1 *input* array of dimension (iter1,nrep1) c rd1 *input* array of dimension (iter1,nrep1) c c____ print the samplers ____________________________________________* parameter (npqmx=10,nmx=400) parameter (nmx2=42230) parameter (iterx=9000) parameter (nrepx=10) parameter (nloopx=20) parameter (nw=500) implicit real*8 (a-h,o-z) common/nconst/np,nq,ixd,npq,npqd,n,nnq,iter1,nrep1,nloop real*8 mphi1(iterx,nrepx,npqmx) real*8 mthe1(iterx,nrepx,npqmx) real*8 mean1(iterx,nrepx) real*8 rd1(iterx,nrepx) real*8 y01(iterx,nrepx,1-npqmx:0) real*8 eps01(iterx,nrepx,1-npqmx:0) real*8 sigma1(iterx,nrepx) c istart=0 iterr=iter1-istart c if (np.ne.0) then do 100 k=1,np sum1=0.0d0 sum2=0.0d0 do 8200 j=1,nrep1 do 8200 i=1+istart,iter1 sum1=sum1+mphi1(i,j,k) 8200 sum2=sum2+mphi1(i,j,k)**2 xmean=sum1/(iterr*nrep1) xvar=(sum2-(iterr*nrep1)*xmean**2)/(iterr*nrep1-1) write (*,*)'mean and std of phi' 100 write (*,9000)xmean,dsqrt(xvar) endif c if (nq.ne.0) then do 200 k=1,nq sum1=0.0d0 sum2=0.0d0 do 8220 j=1,nrep1 do 8220 i=1+istart,iter1 sum1=sum1+mthe1(i,j,k) 8220 sum2=sum2+mthe1(i,j,k)**2 xmean=sum1/(iterr*nrep1) xvar=(sum2-(iterr*nrep1)*xmean**2)/(iterr*nrep1-1) write (*,*)'mean and std of theta' 200 write (*,9000)xmean,dsqrt(xvar) endif c sum1=0.0d0 sum2=0.0d0 do 8230 j=1,nrep1 do 8230 i=1+istart,iter1 sum1=sum1+rd1(i,j) 8230 sum2=sum2+rd1(i,j)**2 xmean=sum1/(iterr*nrep1) xvar=(sum2-(iterr*nrep1)*xmean**2)/(iterr*nrep1-1) write (*,*)'mean and std of d' write (*,9000)xmean,dsqrt(xvar) c sum1=0.0d0 sum2=0.0d0 do 8240 j=1,nrep1 do 8240 i=1+istart,iter1 sum1=sum1+mean1(i,j) 8240 sum2=sum2+mean1(i,j)**2 xmean=sum1/(iterr*nrep1) xvar=(sum2-(iterr*nrep1)*xmean**2)/(iterr*nrep1-1) write (*,*)'mean and std of mean' write (*,9000)xmean,dsqrt(xvar) c sum1=0.0d0 sum2=0.0d0 do 8250 j=1,nrep1 do 8250 i=1+istart,iter1 sum1=sum1+sigma1(i,j) 8250 sum2=sum2+sigma1(i,j)**2 xmean=sum1/(iterr*nrep1) xvar=(sum2-(iterr*nrep1)*xmean**2)/(iterr*nrep1-1) write (*,*)'mean and std of sigma' write (*,9000)xmean,dsqrt(xvar) c if (np.ne.0) then do 110 k=1,np sum1=0.0d0 sum2=0.0d0 do 8260 j=1,nrep1 do 8260 i=1+istart,iter1 sum1=sum1+y01(i,j,k-np) 8260 sum2=sum2+y01(i,j,k-np)**2 xmean=sum1/float(iterr*nrep1) xvar=(sum2-float(iterr*nrep1)*xmean**2)/float(iterr*nrep1-1) write (*,*)'mean and std of y01' 110 write (*,9000)xmean,dsqrt(xvar) endif c if (nq.ne.0) then do 210 k=1,nq sum1=0.0d0 sum2=0.0d0 do 8270 j=1,nrep1 do 8270 i=1+istart,iter1 sum1=sum1+eps01(i,j,k-nq) 8270 sum2=sum2+eps01(i,j,k-nq)**2 xmean=sum1/float(iterr*nrep1) xvar=(sum2-float(iterr*nrep1)*xmean**2)/float(iterr*nrep1-1) write (*,*)'mean and std of eps01' 210 write (*,9000)xmean,dsqrt(xvar) endif c c do 3030 jj=1,nrep1 c do 3030 ii=1+istart,iter1 c030 write (6,1000)(mphi1(ii,jj,k),k=1,np),(mthe1(ii,jj,l),l=1,nq) c030 write (6,1000)(mphi1(ii,jj,l),l=1,np) c030 write (6,1000)(mthe1(ii,jj,l),l=1,nq) c return if (np.ne.0) then do 102 k=1,np c write (6,*)'phi=',k do 102 i=1,iter1 102 write (99,1000)(mphi1(i,j,k),j=1,nrep1) endif c if (nq.ne.0) then do 202 k=1,nq c write (6,*)'theta=',k do 202 i=1,iter1 202 write (99,1000)(mthe1(i,j,k),j=1,nrep1) endif c c write (6,*)'d' do 119 i=1,iter1 119 write (99,1000)(rd1(i,j),j=1,nrep1) c c write (6,*)'mean' do 120 i=1,iter1 120 write (99,1000)(mean1(i,j),j=1,nrep1) c return c write (6,*)'sigma' do 130 i=1,iter1 130 write (99,1000)(sigma1(i,j),j=1,nrep1) c if (np.ne.0) then do 111 k=1,np c write (6,*)'y=',k-np do 111 i=1,iter1 111 write (99,1000)(y01(i,j,k-np),j=1,nrep1) endif c if (nq.ne.0) then do 211 k=1,nq c write (6,*)'eps=',k-nq do 211 i=1,iter1 211 write (99,1000)(eps01(i,j,k-nq),j=1,nrep1) endif c 1000 format(9f8.4) 9000 format(3f10.5) 3033 return end subroutine dblepr return end subroutine intpr return end