* Modified Summer 96 to allow for Full AR/MA matrices/general Sigma * USE MSIMSL **************************************************************** * MVARFIMA PROGRAM * * * ***************************************************************** * modified gibbs sampling for mvarfima(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: array of dim. K,K,np of ar parameters * * THETA: array of dim. K,K,nq of ma parameters * * d: vector ofdim K of frac. diff. parameter * * mu: vector mean parameter * * yo: array of dim.k,np of latent variables * * unobserved historical data * * a0: array of dim. k,nq of latent variables * * unobserved fractional noise * * Sigma ***************************************************************** ***************************************************************** * 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 Wishart * ***************************************************************** ***************************************************************** * 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 * * n: number of data used for the model fitting * * nnq: n+nq * * nk number of series * 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 * * ***************************************************************** * routines called and purpose: * ***************************************************************** * mvarmawt: covariance or correlation of a mvarma model * * fcal: compute log of the posterior density * * to which the complete conditional distributions * * are proportional * * g1: compute the sum of squares * * in the function g1 * * g2: compute the sum of squares * determinant * * in the function g2 * * 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 * * 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 * * upeps: generate sample for (eps0,y0), latent variables * * update: generate sample for (PHI,THETA,d) * * updatemu: generate sample for mean mu * * xomega: compute covariance matrix of (y0,a0) * * (similar to newbold, biometrika,1974) * * and its Cholesky decomposition * * updtsig generate samples of Sigma ****************************************************************** ***************************************************************** * SUBROUTINE lm1 * ****************************************************************** parameter (npqmx=1,nmx=132) parameter (nnqmx=132) parameter (nwmx=500) parameter (nkmx=3) parameter (nrsgmx=12) parameter (nrepx=8) parameter (nloopx=20) implicit real*8 (a-h,o-z) c real sec, CPSEC real*8 datain(nmx,nkmx) real*8 zdata(nkmx,nmx),data(nkmx,nmx) real*8 sigma1(1,nrepx,nkmx,nkmx) real*8 msig1l(1,nrepx,nkmx) real*8 msig1(1,nrepx,nkmx) real*8 mrhol(1,nrepx) real*8 mrho(1,nrepx) real*8 runsigl(nkmx),runrho real*8 runsig(nkmx,nkmx) real*8 stdsig(nkmx),stdrho real*8 mphi1(1,nrepx,nkmx,nkmx,npqmx) real*8 mthe1(1,nrepx,nkmx,nkmx,npqmx) real*8 mean1(1,nrepx,nkmx) real*8 rd1(1,nrepx,nkmx) real*8 y01(1,nrepx,nkmx,1-npqmx:0) real*8 eps01(1,nrepx,nkmx,1-npqmx:0) real*8 runphi(nkmx,nkmx,0:npqmx),runthe(nkmx,nkmx,0:npqmx) real*8 runmn(nkmx),rund(nkmx) ,runeps0(nkmx,1-npqmx:0) real*8 runy0(nkmx,1-npqmx:0),runeps(nkmx,1-npqmx:nmx) real*8 rsig(nrsgmx,nrsgmx),sigfac(nkmx,nkmx) real*8 rsigmu(nkmx,nkmx) integer iseed real*8 grphi1(nrepx,nkmx,nkmx,npqmx) real*8 grthe1(nrepx,nkmx,nkmx,npqmx) real*8 grmn1(nrepx,nkmx),grsig1(nrepx,nkmx,nkmx),grrho1(nrepx) real*8 grsigl1(nrepx,nkmx), grd1(nrepx,nkmx) real*8 temp1(nrepx),temp2(nrepx) real*8 gry01(nrepx,nkmx,1-npqmx:0), greps01(nrepx,nkmx,1-npqmx:0) real*8 grphi2(nrepx,nkmx,nkmx,npqmx) real*8 grthe2(nrepx,nkmx,nkmx,npqmx) real*8 grmn2(nrepx,nkmx),grsig2(nrepx,nkmx,nkmx),grrho2(nrepx) real*8 grsigl2(nrepx,nkmx),grd2(nrepx,nkmx) real*8 gry02(nrepx,nkmx,1-npqmx:0), greps02(nrepx,nkmx,1-npqmx:0) real*8 omgch(npqmx*nkmx,npqmx*nkmx) real*8 cova(nkmx,nkmx,0:nwmx),vi(nkmx,nkmx,0:nnqmx-1) real*8 phi(nkmx,nkmx,0:nnqmx-1,0:nnqmx-1) real*8 s(nkmx,nkmx,1-npqmx-nmx:npqmx-1) real*8 covy(nkmx,nkmx,0:npqmx-1) real*8 tsig(nkmx,nkmx),det(2),work(nkmx) integer ipvt(nkmx) real*4 function dtime EXTERNAL dtime real*4 tarray(2) character*9 today character*8 timestr common/nconst/np,nq,nk,ixd,n,iter1,nrep1,nloop,nrsig,nw common /worksp/rwksp common /grconst/itergr,idf ***************************************************************** * Set up parameter for this subroutine * ***************************************************************** c tarray=dtime() c write(6,*) 'Current time: ', tarray(1),tarray(2) c open(8,file='agg.arma1d0.test.inp') c open(9, file='agg.arma1d0.test.out') call DATE(today) call TIME(timestr) write(6,*) 'Current day/time: ',today,' ',timestr iseed=1653 call rnset(iseed) write(6,*) 'Starting seed: ', iseed read(*,*)n read(*,*)np read(*,*)nq read(*,*)ixd read(*,*)nk read(*,*)nw read(*,*)imnfl read(*,*)isigfl read(*,*)iter1 read(*,*)nrep1 read(*,*)nloop read(*,*)itergr read(*,*)ir write(6,*)'n,np,nq,ixd,nk,nw,imnfl,isigfl,iter1,nrep1,nloop,iterg &r,ir: ',n,np,nq,ixd,nk,nw,imnfl,isigfl,iter1,nrep1, &nloop,itergr,ir npq=np+nq npqd=npq+ixd nnq=n+nq c npk=np*nk nqk=nq*nk npqk=npq*nk ***************************************************************** * initialize runphi(0) and runthe(0) to -I(nk,nk) matrices * ***************************************************************** do ik=1,nk do jk=1,nk runphi(ik,jk,0)=0.0d0 enddo enddo do ik=1,nk runphi(ik,ik,0)=-1.0d0 enddo do ik=1,nk do jk=1,nk runthe(ik,jk,0)=0.0d0 enddo enddo do ik=1,nk runthe(ik,ik,0)=-1.0d0 enddo c Initialize all arrays if (np.ne.0) then c do iitr=1,iter1 do irep=1,nrep1 do ik=1,nk do jk=1,nk do j=1,np mphi1(1,irep,ik,jk,j)=0.0d0 enddo enddo enddo enddo c enddo endif c if(nq.ne.0) then c do iitr=1,iter1 do irep=1,nrep1 do ik=1,nk do jk=1,nk do j=1,nq mthe1(1,irep,ik,jk,j)=0.0d0 enddo enddo enddo enddo c enddo endif c c c if (imnfl .ne. 0) then c do iitr=1,iter1 do irep=1,nrep1 do ik=1,nk mean1(1,irep,ik)=0.0d0 enddo enddo c enddo endif c if (ixd .eq. 1) then c do iitr=1,iter1 do irep=1,nrep1 do ik=1,nk rd1(1,irep,ik)=0.0d0 enddo enddo c enddo endif c do iitr=1,iter1 do irep=1,nrep1 do ik=1,nk do jk=1,nk sigma1(1,irep,ik,jk)=0.0d0 enddo enddo enddo c enddo c if (isigfl .eq.1) then c do iitr=1,iter1 do irep=1,nrep1 mrhol(1,irep)=0.0d0 enddo c enddo c c do iitr=1,iter1 do irep=1,nrep1 do ik=1,nk msig1l(1,irep,ik)=0.0d0 enddo enddo c enddo endif if(np.ne.0) then c do iitr=1,iter1 do irep=1,nrep1 do ik=1,nk do j=1,np y01(1,irep,ik,j-np)=0.0d0 enddo enddo enddo c enddo endif c if (nq.ne.0) then c do iitr=1,iter1 do irep=1,nrep1 do ik=1,nk do j=1,nq eps01(1,irep,ik,j-nq)=0.0d0 enddo enddo enddo c enddo endif ***************************************************************** * * * input initial estimates of parameters and data from a file * * * ***************************************************************** c input initial estimates of mphi(nk,nk,ip), ip=1,....,np * ***************************************************************** if(np.ne.0) then write(6,*)'inimphi' do ip = 1,np do ik=1,nk read(*,*)(mphi1(1,1,ik,jk,ip),jk=1,nk) write(6,1000)(mphi1(1,1,ik,jk,ip),jk=1,nk) enddo enddo endif 1000 format(1x,9f10.6) ***************************************************************** c input initial estimates of theta(nk,nk,iq), iq=1,....,nq ***************************************************************** if (nq.ne.0) then write(6,*)'inimthe' do iq = 1,nq do ik=1,nk read(*,*)(mthe1(1,1,ik,jk,iq),jk=1,nk) write(6,1000)(mthe1(1,1,ik,jk,iq),jk=1,nk) enddo enddo endif ***************************************************************** c input initial estimates of rd1(nk) ***************************************************************** if (ixd.eq.1) then read(*,*) (rd1(1,1,ik),ik=1,nk) write(6,*) 'inird1' write(6,1000)(rd1(1,1,ik),ik=1,nk) endif c set initial estimates of y0 c if(np.ne.0) then do ip=1,np do ik=1,nk y01(1,1,ik,ip-np)=0.0d0 enddo enddo endif c if (nq.ne.0) then do ik=1,nk do iq=1,nq eps01(1,1,ik,iq-nq)=0.0d0 enddo enddo endif ***************************************************************** c input initial estimates of mean1(nk) ***************************************************************** if (imnfl .ne.0) then write(6,*)'inimean' read(*,*)(mean1(1,1,ik),ik=1,nk) write(6,1000)(mean1(1,1,ik),ik=1,nk) else do ik=1,nk mean1(1,1,ik)=0.0d0 enddo endif ***************************************************************** * input initial estimates of sigma :cov. matrix of errors * ***************************************************************** write(6,*)'inisigma' if (isigfl.eq.0) then do ik=1,nk read(*,*)(sigma1(1,1,ik,jk),jk=1,nk) write(6,1000)(sigma1(1,1,ik,jk),jk=1,nk) enddo do ik=1,nk do jk=1,nk runsig(ik,jk)=sigma1(1,1,ik,jk) enddo enddo c Find Factorization of sigma1 (needed to generate samples from a Wishart) call DCHFAC(nk,runsig,nkmx,100.0d0*DMACH(4),irank,sigfac,nkmx) else ***************************************************************** * for concurrent case, input res. variances for nk series, * * transform them * ***************************************************************** do ik=1,nk read(*,*)sigma1(1,1,ik,ik) msig1(1,1,ik)=sigma1(1,1,ik,ik) msig1l(1,1,ik)=dlog(dsqrt(sigma1(1,1,ik,ik))) write(6,1001)sigma1(1,1,ik,ik),msig1l(1,1,ik) enddo 1001 format(1x,2f10.4) ***************************************************************** * for concurrent case, input initial estimate of rho * * transform by Fisher-type transformation * ***************************************************************** read(*,*) rhoini mrho(1,1)=rhoini mrhol(1,1)=dlog(rhoini/(1.0-rhoini)) write(6,1002)'inirho,rhofish',rhoini,mrhol(1,1) 1002 format(1x,a,2f10.4) ***************************************************************** * for concurrent case, input initial estimate of stddev of * * series variances * do ik=1,nk read(*,*) stdsig(ik) write(6,1003)'stdsig',stdsig(ik) enddo 1003 format(1x,a,f10.4) * Read Stdrho read (*,*) stdrho write(6,1003)'stdrho',stdrho * form the cov. matrix do ik=1,nk do jk=1,nk runsig(ik,jk)=sigma1(1,1,ik,ik)*sigma1(1,1,jk,jk) if(ik.ne.jk) runsig(ik,jk)=runsig(ik,jk)*rhoini c write(6,*)'runsig ',runsig(ik,jk) enddo enddo endif ***************************************************************** ***************************************************************** * input initial estimates of covariance between Phi Theta D(Upper Triangle) * into rsig and convert to its Cholesky decomposition * ****************************************************************** nrsig=npq*(nk**2)+nk*ixd do ik = 1, nrsig read(*,*)( rsig(ik,jk),jk=1,nrsig) c write(6,*) (rsig(ik,jk),jk=1,nrsig) enddo call DCHFAC(nrsig,rsig,nrsgmx,100*DMACH(4),irank,rsig,nrsgmx) ***************************************************************** * input initial estimates of covariance of mu(nk) (Upper Triangle) * into rsigmu and convert rsigmu to its cholesky decomposition* ***************************************************************** if (imnfl .ne. 0) then do ik=1,nk read(*,*) (rsigmu(ik,jk),jk=1,nk) enddo call DCHFAC(nk,rsigmu,nkmx,100*DMACH(4),irank,rsigmu,nkmx) endif c set initial estimates of y0 and eps0 c if(np.ne.0) then do ip=1,np do ik=1,nk y01(1,1,ik,ip-np)=0.0d0 enddo enddo endif c if (nq.ne.0) then do ik=1,nk do iq=1,nq eps01(1,1,ik,iq-nq)=0.0d0 enddo enddo endif ***************************************************************** * input data in datain and read into zdata * ***************************************************************** do i=1,n read(*,*)(datain(i,ik),ik=1,nk) enddo do ik=1,nk do i=1,n zdata(ik,i)=datain(i,ik) enddo enddo if(np.ne.0) then do ip=1,np do ik=1,nk do jk=1,nk runphi(ik,jk,ip)=mphi1(1,1,ik,jk,ip) enddo enddo enddo endif if (nq.ne.0) then do iq=1,nq do ik=1,nk do jk=1,nk runthe(ik,jk,iq)=mthe1(1,1,ik,jk,iq) enddo enddo enddo endif c if(ixd.eq.1) then do ik=1,nk rund(ik)=rd1(1,1,ik) enddo endif c c obtain cholesky decomposition of omega, put into omgch if(npq.ne.0) then call xomega(omgch,runphi,runthe,rund,runsig) c Find Factorization of omega (needed to generate samples of y0,e0) call DCHFAC(nk*(np+nq),omgch,nkmx*npqmx,100.0d0*DMACH(4), & irank,omgch,nkmx*npqmx) endif c ***************************************************************** * perturb the initial values to obtain * * starting values for nrep1 parallel gibbs chains. * ***************************************************************** if (nrep1 .ge. 1) / call perturb (mphi1,mthe1,y01,eps01,sigma1,mean1,rd1, / rsig,rsigmu,omgch,sigfac,msig1l,mrhol, / stdsig,stdrho,imnfl,isigfl) ********************************************************************** ***************************************************************** * begin nrep1 independent Gibbs chains * ***************************************************************** do 1120 jjrep=1,nrep1 write(6,*) 'Replication #:', jjrep c sec=dble(CPSEC()) ***************************************************************** c Initialize sums for computation of Gelman-Rubin ****************************************************************** c if (np .ne. 0) then do k=1,np do ik=1,nk gry01(jjrep,ik,k-np)=0.0d0 gry02(jjrep,ik,k-np)=0.0d0 do jk=1,nk grphi1(jjrep,ik,jk,k)=0.0d0 grphi2(jjrep,ik,jk,k)=0.0d0 enddo enddo enddo endif c if (nq .ne. 0) then do k=1,nq do ik=1,nk greps01(jjrep,ik,k-nq)=0.0d0 greps02(jjrep,ik,k-nq)=0.0d0 do jk=1,nk grthe1(jjrep,ik,jk,k)=0.0d0 grthe2(jjrep,ik,jk,k)=0.0d0 enddo enddo enddo endif c if (ixd .eq. 1) then do ik=1,nk grd1(jjrep,ik)=0.0d0 grd2(jjrep,ik)=0.0d0 enddo endif c if (imnfl .eq. 1) then do ik=1,nk grmn1(jjrep,ik)=0.0d0 grmn2(jjrep,ik)=0.0d0 enddo endif c if(isigfl .eq.0) then do ik=1,nk do jk=1,nk grsig1(jjrep,ik,jk)=0.0d0 grsig2(jjrep,ik,jk)=0.0d0 enddo enddo endif c if(isigfl.eq.1) then do ik=1,nk grsigl1(jjrep,ik)=0.0d0 grsigl2(jjrep,ik)=0.0d0 enddo c grrho1(jjrep)=0.0d0 grrho2(jjrep)=0.0d0 endif ***************************************************************** * set starting value of each chain * ***************************************************************** if(nq.ne.0) then do j=1,nq do ik=1,nk runeps0(ik,j-nq)=eps01(1,jjrep,ik,j-nq) write(6,*)'runeps0',runeps0(ik,j-nq) enddo enddo endif c if(np.ne.0) then do j=1,np do ik=1,nk runy0(ik,j-np)= y01(1,jjrep,ik,j-np) write(6,*)'runy0',runy0(ik,j-np) enddo enddo endif if(np.ne.0) then do j=1,np do ik=1,nk do jk=1,nk runphi(ik,jk,j)=mphi1(1,jjrep,ik,jk,j) enddo enddo enddo do j=1,np do ik=1,nk runy0(ik,1-j)=0.0d0 enddo enddo endif c write runphi do j=1,np write(6,95)((runphi(ik,jk,j), ik=1,nk),jk=1,nk) enddo if (nq.ne.0) then do j=1,nq do ik=1,nk do jk=1,nk runthe(ik,jk,j)=mthe1(1,jjrep,ik,jk,j) enddo enddo enddo do j=1,nq do ik=1,nk runeps0(ik,1-j)=0.0d0 enddo enddo endif c write runthe do j=1,nq write(6,95)((runthe(ik,jk,j), ik=1,nk),jk=1,nk) enddo c write rund if (ixd.eq.1) then do ik=1,nk rund(ik)=rd1(1,jjrep,ik) write(6,*)'rund',rund(ik) enddo else do i=1,nk rund(ik)=0.0d0 enddo endif c *********************************************************************** * Setup Sigma if (isigfl.eq.0) then do ik=1,nk do jk=1,nk runsig(ik,jk)=sigma1(1,jjrep,ik,jk) write(6,1003)'runsig ',runsig(ik,jk) enddo enddo else do ik=1,nk runsigl(ik)=msig1l(1,jjrep,ik) enddo write(6,1000)(runsigl(ik),ik=1,nk) runrho=mrhol(1,jjrep) write(6,1003)'runrho',runrho * form the untransformed cov. matrix each time rho=dexp(runrho)/(1.0d0+dexp(runrho)) do ik=1,nk do jk=1,nk runsig(ik,jk)=dexp(runsigl(ik))*dexp(runsigl(jk)) if(ik.ne.jk) runsig(ik,jk)=runsig(ik,jk)*rho enddo enddo endif tsig=runsig call DGEFA(tsig,nkmx,nk,ipvt,info) call DGEDI(tsig,nkmx,nk,ipvt,det,work,10) desig=dlog(det(1)*10.0d0**det(2)) if (imnfl .ne.0) then do ik=1,nk runmn(ik)=mean1(1,jjrep,ik) write(6,1003)'runmn',runmn(ik) enddo else do ik=1,nk runmn(ik)=0.0d0 enddo endif ***************************************************************** * subtracts mean from zdata to form data * ***************************************************************** do ik=1,nk do in=1,n data(ik,in)=zdata(ik,in)-runmn(ik) enddo enddo c **************************************************************** * compute a* vector = (a_0, y_0) * ***************************************************************** call musig(runphi,runthe,runy0,runeps0,data,runeps) ***************************************************************** * compute f first time * ************************************************************ call gcova (cova,rund,runsig) call gc(cova,det1,vi,phi) call g1 (runeps,phi,vi,a1) c call g1test(cova,runeps,det1t,a1t) if(np.ne.0) then call gsy(cova,runphi,runthe,s,covy) call g2(phi,vi,runeps,runy0,covy,s,a2,de2) endif call fcal(f,a1,a2,det1,de2,desig) c write(6,*) 'Fcal first time' c write(6,*) 'f,a1,a2,det1,de2,desig: ', f,a1,a2,det1,de2,desig c call fcal(f,a1t,a2,det1t,de2,desig) c write(6,*) 'f,a1t,det1t,desig: ', f,a1t,det1t,desig c c sec=CPSEC() c write(6,*) 'Iteration 0: CPU secs: ', sec ***************************************************************** * run MCMC with iter1-1 iterations * ***************************************************************** do 1100 iitr=1,iter1-1 write (6,*)'MC Iteration: ', iitr c * wall_seconds()-time_now ***************************************************************** * subtracts mean from zdata to form data * ***************************************************************** do ik=1,nk do in=1,n data(ik,in)=zdata(ik,in)-runmn(ik) enddo enddo c ***************************************************************** * update (PHI,THETA,d) * ***************************************************************** call update (f,phi,vi,runphi,runthe,runy0,runeps0, / runeps,runsig,rsig,data,rund,cova,covy,s, / a1,a2,det1,de2,desig) c write(6,*) 'after update:f,a1,a2,det1,de2,desig: ', c & f,a1,a2,det1,de2,desig c c write current output (every r th after itergr iterations) c for phi if ((np.ne.0) .and. (iitr .ge. itergr)) then if (mod(iitr,ir).eq.0) then do j=1,np write(6,95)((runphi(ik,jk,j),ik=1,nk),jk=1,nk) enddo 95 format(9f8.4) endif c update running totals after itergr iterations do ik=1,nk do jk=1,nk do k=1,np grphi1(jjrep,ik,jk,k)=grphi1(jjrep,ik,jk,k)+ & runphi(ik,jk,k) grphi2(jjrep,ik,jk,k)=grphi2(jjrep,ik,jk,k)+ & runphi(ik,jk,k)**2 enddo enddo enddo endif c c write current output (every r th after itergr iterations) c for theta if ((nq.ne.0) .and. (iitr .ge. itergr)) then if (mod(iitr,ir).eq.0) then do j=1,nq write(6,95)((runthe(ik,jk,j),ik=1,nk),jk=1,nk) enddo endif c update running totals after itergr iterations do ik=1,nk do jk=1,nk do k=1,nq grthe1(jjrep,ik,jk,k)=grthe1(jjrep,ik,jk,k)+ & runthe(ik,jk,k) grthe2(jjrep,ik,jk,k)=grthe2(jjrep,ik,jk,k)+ & runthe(ik,jk,k)**2 enddo enddo enddo endif c c write current output (every r th after itergr iterations) c for d if ((ixd.ne.0) .and. (iitr .ge. itergr)) then if (mod(iitr,ir).eq.0) then write(6,95)(rund(ik), ik=1,nk) endif c update running totals after itergr iterations do ik=1,nk grd1(jjrep,ik)=grd1(jjrep,ik)+ & rund(ik) grd2(jjrep,ik)=grd2(jjrep,ik)+ & rund(ik)**2 enddo endif ***************************************************************** * update latent variables * ***************************************************************** if(npq.ne.0) then call upeps (f,phi,vi,runphi,runthe,runy0,runeps0, 1 omgch,data,runeps,covy,s,a1,a2, 1 det1,de2,desig) c write(6,*) 'after upeps:f,a1,a2,det1,de2,desig: ', c & f,a1,a2,det1,de2,desig c write current output and update running totals if ((np.ne.0) .and. (iitr .ge. itergr)) then if (mod(iitr,ir).eq.0) then do k=1,np write(6,1000)(runy0(ik,k-np),ik=1,nk) enddo endif do k=1,np do ik=1,nk gry01(jjrep,ik,k-np)=gry01(jjrep,ik,k-np)+ & runy0(ik,k-np) gry02(jjrep,ik,k-np)=gry02(jjrep,ik,k-np)+ & runy0(ik,k-np)**2 enddo enddo endif c if ((nq.ne.0) .and. (iitr .ge. itergr)) then if (mod(iitr,ir).eq.0) then do k=1,nq write(6,95)(runeps0(ik,k-nq),ik=1,nk) enddo endif do k=1,nq do ik=1,nk greps01(jjrep,ik,k-nq)=greps01(jjrep,ik,k-nq)+ & runeps0(ik,k-nq) greps02(jjrep,ik,k-nq)=greps02(jjrep,ik,k-nq)+ & runeps0(ik,k-nq)**2 enddo enddo endif endif c ***************************************************************** * update mu * ***************************************************************** if (imnfl .ne. 0) then call updatemu (f,phi,vi,runmn,runphi,runthe,runy0,runeps0, 1 runeps,rsigmu,data,covy,s,a1,a2, 1 det1,de2,desig) c write(6,*) 'after updatemu:f,a1,a2,det1,de2,desig: ', c & f,a1,a2,det1,de2,desig c write current output and update running totals if (iitr. ge. itergr) then if(mod(iitr,ir).eq.0) then write(6,95)(runmn(ik),ik=1,nk) endif c do ik=1,nk grmn1(jjrep,ik)=grmn1(jjrep,ik)+runmn(ik) grmn2(jjrep,ik)=grmn2(jjrep,ik)+runmn(ik)**2 enddo endif endif ***************************************************************** * update sigma * ***************************************************************** if (isigfl .eq.0) then call updtsig(f,phi,vi,runphi,runthe,runy0, 1 runeps,runsig,sigfac,rund,cova,covy,s,a1,a2, 1 det1,de2,desig) c write(6,*) 'after updtsig:f,a1,a2,det1,de2,desig: ', c & f,a1,a2,det1,de2,desig c write current output and running totals if (iitr .ge.itergr) then if(mod(iitr,ir).eq.0) then do ik =1,nk write(6,1000)(runsig(ik,jk),jk=1,nk) enddo endif c do ik=1,nk do jk=1,nk grsig1(jjrep,ik,jk)=grsig1(jjrep,ik,jk)+runsig(ik,jk) grsig2(jjrep,ik,jk)=grsig2(jjrep,ik,jk)+runsig(ik,jk)**2 enddo enddo endif ***************************************************************** * update sigma and rho * ***************************************************************** c write(6,*) 'calling updtsig' else call updtsig2(f,phi,vi,runphi,runthe,stdsig,runy0, 1 runeps,runsigl,rund,runrho,cova,covy,s,a1,a2, 1 det1,de2,desig) c write current output and running totals if (iitr .ge. itergr) then if(mod(iitr,ir).eq.0) then write(6,1000)(dexp(runsigl(ik)),ik=1,nk) c write(6,1000)(runsigl(ik),ik=1,nk) endif c do ik=1,nk grsigl1(jjrep,ik)=grsigl1(jjrep,ik)+(dexp(runsigl(ik)))**2 grsigl2(jjrep,ik)=grsigl2(jjrep,ik)+((dexp(runsigl(ik)))**2)**2 enddo endif c write(6,*)'after updtsig f,a1,desig', f,a1,desig c call updtrho(f,phi,vi,runphi,runthe,stdrho,runy0, 1 runeps,runsigl,runrho ,rund,cova,covy,s,a1,a2, 1 det1,de2,desig) c write current output and running totals rho=dexp(runrho)/(1.0d0+dexp(runrho)) if (iitr .ge. itergr) then if(mod(iitr,ir).eq.0) then write(6,9850) rho endif grrho1(jjrep)=grrho1(jjrep)+rho grrho2(jjrep)=grrho2(jjrep)+rho**2 endif c c write(6,*)'after updtrho f,a1,desig', f,a1,desig c compute runsig do ik=1,nk do jk=1,nk runsig(ik,jk)=(dexp(runsigl(ik)))*(dexp(runsigl(jk))) if(ik.ne.jk) runsig(ik,jk)=runsig(ik,jk)*rho enddo enddo endif ***************************************************************** ****************************************************************** c sec=CPSEC() c write(6,*) 'Iteration ', iitr, ' CPU secs: ', sec 1100 continue 1120 continue ***************************************************************** ********************************************************************** * Compute G-R statistics ********************************************************************** if (nrep1 .gt. 1) then c for each phi parameter if (np .ne. 0) then do k=1,np do ik=1,nk do jk=1,nk do jjrep=1,nrep1 temp1(jjrep)=grphi1(jjrep,ik,jk,k) temp2(jjrep)=grphi2(jjrep,ik,jk,k) enddo call grubin(temp1,temp2,gr) write(6,*)'ip,gr for phi',k,gr enddo enddo enddo endif c for each theta parameter if (nq .ne. 0) then do k=1,nq do ik=1,nk do jk=1,nk do jjrep=1,nrep1 temp1(jjrep)=grthe1(jjrep,ik,jk,k) temp2(jjrep)=grthe2(jjrep,ik,jk,k) enddo call grubin(temp1,temp2,gr) write(6,*)'iq,gr for theta',k,gr enddo enddo enddo endif c for d parameters if (ixd.ne. 0) then do ik=1,nk do jjrep=1,nrep1 temp1(jjrep)=grd1(jjrep,ik) temp2(jjrep)=grd2(jjrep,ik) enddo call grubin(temp1,temp2,gr) write(6,*)'gr for d ', ik, gr enddo endif c for mean parameters if (imnfl.ne. 0) then do ik=1,nk do jjrep=1,nrep1 temp1(jjrep)=grmn1(jjrep,ik) temp2(jjrep)=grmn2(jjrep,ik) enddo call grubin(temp1,temp2,gr) write(6,*)'gr for mean ', ik, gr enddo endif if (isigfl.eq.0) then c for general sigma matrix do ik=1,nk do jk=1,nk do jjrep=1,nrep1 temp1(jjrep)=grsig1(jjrep,ik,jk) temp2(jjrep)=grsig2(jjrep,ik,jk) enddo call grubin(temp1,temp2,gr) write(6,*)'gr for sigma',ik,jk,gr enddo enddo else c for sigma parameters and rho do ik=1,nk do jjrep=1,nrep1 temp1(jjrep)=grsigl1(jjrep,ik) temp2(jjrep)=grsigl2(jjrep,ik) c write(6,*)'temp1 for sigma',temp1(jjrep) c write(6,*)'temp2 for sigma',temp2(jjrep) enddo call grubin(temp1,temp2,gr) write(6,*)'gr for sigma', ik,gr enddo c do jjrep=1,nrep1 temp1(jjrep)=grrho1(jjrep) temp2(jjrep)=grrho2(jjrep) enddo call grubin(temp1,temp2,gr) write(6,*)'gr for rho', gr endif c endif c write output call prx0 (grphi1,grphi2,grthe1,grthe2,grd1,grd2,gry01,gry02, &greps01,greps02,grmn1,grmn2,grsig1,grsig2,grsigl1,grsigl2, &grrho1,grrho2,imnfl,isigfl) c write random seed call rnget(iseed) write(6,*)'seed',iseed 9700 format(1x,3f10.4) 9850 format(1x,f10.4) c tarray=dtime() c write(*,*) 'Current time: ', tarray(1),tarray(2) call DATE(today) call TIME(timestr) write(6,*) 'Current day/time: ',today,' ',timestr STOP END ********************************************************************** * * * SUBROUTINE fcal * * * ********************************************************************** subroutine fcal (f,a1,a2,det1,de2,desig) ********************************************************************** * 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 * * det1 *input * scalar. determinant in the g1 function * * de2 *input * scalar. determinant in the g2 function * * detsig*input * scalar. contains determinant of WN cov. SIGMA * implicit real*8 (a-h,o-z) common/nconst/np,nq,nk,ixd,n,iter1,nrep1,nloop,nrsig,nw c if(np.eq.0)then f=-0.5d0*(det1+a1+desig) else f=-0.5d0*(det1+de2+a1+a2+desig) endif return end ********************************************************************** * * * SUBROUTINE g1 * * * ********************************************************************** subroutine g1 (runeps,phi,v,a1) ********************************************************************** c use: compute sum of squares term in g1 function c parameters of routine: c runeps *input * array of size (I,1-nq:n).-- Please note: the order of c elements in runeps is (eps_1,{1-nq},...,eps_I,{1-nq},.....,eps_1,n,... c eps_I,n c phi *input * array of dimension (I,I,0:n+nq-1,0:n+nq-1) contains c forward partial correlation matrices needed for c sigma 22 inverse c v *input * array of dimension (I,I,0:n+nq-1) contains block c diagonal elements needed for sigma 22 inverse c a1 *output* scalar. computed value of sum of squares term c in the g1 function c c other Routines Used: MATMUL,MATTRAN,VMATMUL parameter (npqmx=1,nmx=132) parameter (nnqmx=132) parameter (nkmx=3) implicit real*8 (a-h,o-z) integer i,j,k,l,m real*8 runeps(nkmx,1-npqmx:nmx),temp1(nkmx,0:nnqmx-1) real*8 phi(nkmx,nkmx,0:nnqmx-1,0:nnqmx-1) real*8 v(nkmx,nkmx,0:nnqmx-1) common/nconst/np,nq,nk,ixd,n,iter1,nrep1,nloop,nrsig,nw c c Multiply runeps and beta do i=0,n+nq-1 do l=1,nk c sum=sum+phi1(i-j,i)*runeps(1-nq+j) temp1(l,i)=0.0d0 do j=0,i do m=1,nk temp1(l,i)=temp1(l,i)+phi(l,m,i-j,i)*runeps(m,1-nq+j) enddo enddo enddo enddo c Multiply temp1' with V with temp1 a1=0.0d0 do i=0,n+nq-1 c sum=sum+temp1(i)*V(i)*temp1(i) do j=1,nk do k=1,nk a1=a1+temp1(j,i)*v(j,k,i)*temp1(k,i) enddo enddo enddo return end ********************************************************************* subroutine g2(phi,v,runeps,runy0,covy,s,a2,de2) ********************************************************************* c This subroutine computes the LN of determinant and sum of squares in the c g2 function, i.e. sigma_11-sigma_12 sigma_22 ^-1 sigma_2,1 c parameters: c phi :input array of size (nkmx,nkmx,0:n+nq-1,0:n+nq-1) Contains Forward c v :input array of size (nkmx,nkmx,0:n+nq-1) Contains v_j^-1 c runeps:input array of size (nkmx,1-nq:n) c runy0 :input array of size (nkmx,nkmx,1-npqmx:0) c covy :input array of size (nkmx,nkmx,0:npqmx-1) Covariance of ARFIMA(p,d, c s :input array of size (nkmx,nkmx,1-np-n:nq-1) Contains eta elements c a2 :output scaler containing sum-of-squares value c de2 :output scaler containing LN of determinant c Other routines used: MATMUL,MATTRAN,VMATMUL,SCALER,DGEFA,DGEDI parameter (npqmx=1,nmx=132) parameter (nnqmx=132) parameter (nkmx=3) implicit real*8 (a-h,o-z) integer i,j,k,l,m,o real*8 runy0(nkmx,1-npqmx:0) real*8 runeps(nkmx,1-npqmx:nmx) real*8 v(nkmx,nkmx,0:nnqmx-1) real*8 phi(nkmx,nkmx,0:nnqmx-1,0:nnqmx-1) real*8 s(nkmx,nkmx,1-npqmx-nmx:npqmx-1) real*8 covy(nkmx,nkmx,0:npqmx-1) real*8 a2,de2,det(2),work(npqmx*nkmx) integer ipvt(npqmx*nkmx) double precision temp1(nkmx,nkmx,0:nnqmx-1,npqmx),w1(nkmx,nkmx), & w2(nkmx,nkmx),w3(nkmx,nkmx),w4(nkmx,nkmx),w5(nkmx),w6(nkmx), & sum(nkmx,nkmx),zero(nkmx,nkmx),sumv(nkmx),zerov(nkmx), & temp2(npqmx*nkmx,npqmx*nkmx),vmu(nkmx*npqmx), & temp3(nkmx,0:nnqmx-1),fac(npqmx*nkmx,npqmx*nkmx) common/nconst/np,nq,nk,ixd,n,iter1,nrep1,nloop,nrsig,nw c write(6,*) 'in g2' call scaler(nk,0.0d0,zero) do i=1,nk zerov(i)=0.0d0 enddo c Multiply Beta and sigma_21' do i=0,n+nq-1 do j=1,np c Initialize temp1 do l=1,nk do m=1,nk temp1(l,m,i,j)=0.0d0 enddo enddo do k=0,i c sum=sum+phi(i-k,i)*s((j-i)-(np-nq))' do l=1,nk do m=1,nk do o=1,nk temp1(l,m,i,j)=temp1(l,m,i,j)+ & phi(l,o,i-k,i)*s(m,o,(j-k-1)-(np-nq)) enddo enddo enddo enddo enddo enddo c Multiply temp1' with V with temp1 c and subtract sigma_12 *sigma_22 inv * sigma_21 from sigma_11 do i=1,np do j=1,np do ik=1,nk do jk=1,nk sum(ik,jk)=zero(ik,jk) enddo enddo do k=0,n+nq-1 c Multiplication of I*I block do l=1,nk do m=1,nk w1(l,m)=temp1(m,l,k,i) w2(l,m)=v(l,m,k) w3(l,m)=temp1(l,m,k,j) enddo enddo call MATMUL(w1,w2,w4,nk,nk,nk) call MATMUL(w4,w3,w1,nk,nk,nk) do ik=1,nk do jk=1,nk sum(ik,jk)=sum(ik,jk)+w1(ik,jk) enddo enddo c sum=sum+w1 enddo if ((i-j) .ge.0) then do l=1,nk do m=1,nk temp2((i-1)*nk+l,(j-1)*nk+m)=covy(l,m,(i-j))-sum(l,m) enddo enddo else do l=1,nk do m=1,nk temp2((i-1)*nk+l,(j-1)*nk+m)=covy(m,l,(j-i))-sum(l,m) enddo enddo endif enddo enddo c Find determinant and inverse of temp2 c write(6,*) 'temp2' c do i=1,nk*np c write(6,*) (temp2(i,j),j=1,np*nk) c enddo call dgefa(temp2,npqmx*nkmx,np*nk,ipvt,info) call dgedi(temp2,npqmx*nkmx,np*nk,ipvt,det,work,11) c call DLFTDS(np*nk,temp2,npqmx*nkmx,fac,npqmx*nkmx) c call DLFDDS(np*nk,fac,npqmx*nkmx,det1,det2) c call DLINDS(np*nk,temp2,npqmx*nkmx,temp2,npqmx*nkmx) c write(6,*) 'det: ', det1*(10.**det2) de2=dlog(det(1)*(10.0d0**det(2)) ) c To Compute vmu -- Multiply sigma_12 with sigma_22 inv with runeps c First multiply V and Beta and runeps do i=0,n+nq-1 do ik=1,nk sumv(ik)=zerov(ik) enddo c sumv=zerov do j=0,i do l=1,nk do m=1,nk w2(l,m)=phi(l,m,i-j,i) w1(l,m)=v(l,m,i) enddo w5(l)=runeps(l,1-nq+i) enddo call MATMUL(w1,w2,w3,nk,nk,nk) call VMATMUL(w3,w5,w6,nk,nk) do ik=1,nk sumv(ik)=sumv(ik)+w6(ik) enddo c sumv=sumv+w6 enddo do l=1,nk temp3(l,i)=sumv(l) enddo enddo c Then Multiply result by sigma_12 and Beta' (temp1') do i=1,np do ik=1,nk sumv(ik)=zerov(ik) enddo c sumv=zerov do j=0,n+nq-1 do l=1,nk do m=1,nk sumv(l)=sumv(l)+temp1(m,l,j,i)*temp3(m,j) enddo enddo enddo do l=1,nk vmu((i-1)*nk+l)=runy0(l,-np+i)-sumv(l) enddo enddo c Now compute SSQ a2=0.0d0 do i=1,nk*np do j=1,nk*np a2=a2+vmu(i)*temp2(i,j)*vmu(j) enddo enddo c c write(6,*)'de2 and a2 in g2',de2,a2 return end ********************************************************************** * * * SUBROUTINE gc * * * ********************************************************************** subroutine gc(cova,det1,vi,phi) ********************************************************************** c Compute determinant of sigma_{2,2} and elements needed in factorization c of inverse of sigma_{2,2} c Input: cova -- array of size(I,I,0:nw) containing covariances of MVARFIMA( c Output:det1 -- scaler containing determinant of sigma_{2,2} c vi -- array of size (I,I,0:n+q-1) containing v_j^{-1} c phi -- array of size (I,I,0:n+nq-1,0:n+nq-1) containing forward c partial correlation matrices (needed for sigma_{2,2}^-1 parameter(nkmx=3,nmx=132,npqmx=1,nwmx=500) parameter(nnqmx=132) implicit real*8 (a-h,o-z) real*8 v(nkmx,nkmx),vinv(nkmx,nkmx),vstrinv(nkmx,nkmx), & vstr(nkmx,nkmx),vi(nkmx,nkmx,0:nnqmx-1) real*8 del(nkmx,nkmx),delstr(nkmx,nkmx) real*8 phi(nkmx,nkmx,0:nnqmx-1,0:nnqmx-1), & phistr(nkmx,nkmx,0:nnqmx-1,0:nnqmx-1) real*8 w1(nkmx,nkmx),w2(nkmx,nkmx),w3(nkmx,nkmx) real*8 w4(nkmx,nkmx),w5(nkmx,nkmx),w6(nkmx,nkmx) real*8 cova(nkmx,nkmx,0:nwmx),det1,det(2),work(nkmx) integer ipvt(nkmx) common/nconst/np,nq,nk,ixd,n,iter1,nrep1,nloop,nrsig,nw c Initialize phi(0,j) to Identity c call scaler(nk,1.0d0,w1) c do j=0,n+nq-1 do l=1,nk do m=1,nk phi(l,m,0,j)=w1(l,m) phistr(l,m,0,j)=w1(l,m) enddo enddo enddo do l=1,nk do m=1,nk v(l,m)=cova(l,m,0) vstr(l,m)=v(l,m) del(l,m)=cova(l,m,1) delstr(l,m)=cova(m,l,1) enddo enddo c Compute inverse of v and determinant of v (inverse returned in v) call DGEFA(v,nkmx,nk,ipvt,info) call DGEDI(v,nkmx,nk,ipvt,det,work,11) det1=dlog(det(1)*(10.0d0**det(2))) c write(6,*) 'in gc: 1st', det1 do l=1,nk do m=1,nk vi(l,m,0)=v(l,m) vinv(l,m)=v(l,m) vstrinv(l,m)=v(l,m) enddo enddo c do j=1,n+nq-1 do ik=1,nk do jk=1,nk w1(ik,jk)=-del(ik,jk) enddo enddo call MATMUL(w1,vstrinv,w2,nk,nk,nk) do ik=1,nk do jk=1,nk w1(ik,jk)=-delstr(ik,jk) enddo enddo c w1=-delstr call MATMUL(w1,vinv,w3,nk,nk,nk) do l=1,nk do m=1,nk phi(l,m,j,j)=w2(l,m) phistr(l,m,j,j)=w3(l,m) enddo enddo do i=1,j-1 do l=1,nk do m=1,nk w1(l,m)=phi(l,m,j,j) w2(l,m)=phistr(l,m,j-i,j-1) w4(l,m)=phistr(l,m,j,j) w5(l,m)=phi(l,m,j-i,j-1) enddo enddo call MATMUL(w1,w2,w3,nk,nk,nk) call MATMUL(w4,w5,w6,nk,nk,nk) do l=1,nk do m=1,nk phi(l,m,i,j)=phi(l,m,i,j-1)+w3(l,m) phistr(l,m,i,j)=phistr(l,m,i,j-1)+w6(l,m) enddo enddo enddo do l=1,nk do m=1,nk v(l,m)=cova(l,m,0) vstr(l,m)=v(l,m) del(l,m)=cova(l,m,j+1) delstr(l,m)=cova(m,l,j+1) enddo enddo do i=1,j do l=1,nk do m=1,nk w1(l,m)=phi(l,m,i,j) w2(l,m)=cova(l,m,i) w5(l,m)=phistr(l,m,i,j) w6(l,m)=cova(l,m,j+1-i) enddo enddo call MATTRAN(w2,w3,nk,nk) call MATMUL(w1,w3,w4,nk,nk,nk) do ik=1,nk do jk=1,nk v(ik,jk)=v(ik,jk)+w4(ik,jk) enddo enddo c v=v+w4 call MATMUL(w5,w2,w3,nk,nk,nk) do ik=1,nk do jk=1,nk vstr(ik,jk)=vstr(ik,jk)+w3(ik,jk) enddo enddo c vstr=vstr+w3 call MATMUL(w1,w6,w3,nk,nk,nk) do ik=1,nk do jk=1,nk del(ik,jk)= del(ik,jk)+w3(ik,jk) enddo enddo c del=del+w3 call MATTRAN(w6,w3,nk,nk) call MATMUL(w5,w3,w4,nk,nk,nk) do ik=1,nk do jk=1,nk delstr(ik,jk)= delstr(ik,jk)+w4(ik,jk) enddo enddo c delstr=delstr+w4 enddo call DGEFA(v,nkmx,nk,ipvt,info) call DGEDI(v,nkmx,nk,ipvt,det,work,11) det1=det1+dlog(det(1)*(10.0d0**det(2))) c write(6,*) 'in gc: ', det(1)*(10.0d0**det(2)) call DGEFA(vstr,nkmx,nk,ipvt,info) call DGEDI(vstr,nkmx,nk,ipvt,det,work,01) do l=1,nk do m=1,nk vi(l,m,j)=v(l,m) vinv(l,m)=v(l,m) vstrinv(l,m)=vstr(l,m) enddo enddo enddo return end ********************************************************************** * * * SUBROUTINE gcova * * * ********************************************************************** subroutine gcova(cova,rd,sigma) ********************************************************************** c Computes the covariance of a MVARFIMA(0,d_i,0) process c See Sowell, 1989. c Input: rd -- vector of length I containing current values of d c sigma -- array of size (IxI) containing var-cov of errors c Output: cova -- array of size (I,I,0:nw) containing acf c Other routines used: GAMA (subroutine to compute GAMMA function) parameter(nkmx=3,nmx=132,npqmx=1,nwmx=500) parameter(nnqmx=132) implicit real*8 (a-h,o-z) real*8 sigma(nkmx,nkmx),rd(nkmx),cova(nkmx,nkmx,0:nwmx) integer nk,nw,i,j,k common/nconst/np,nq,nk,ixd,n,iter1,nrep1,nloop,nrsig,nw c do i=1,nk do j=1,nk cova(i,j,0)=sigma(i,j)*DGAMMA(1.0d0-rd(i)-rd(j))/ & (DGamma(1.0d0-rd(j))*DGamma(1.0d0-rd(i))) enddo enddo do i=1,nk do j=1,nk do k=1,nw cova(i,j,k)=-cova(i,j,k-1)* & (1.d0-rd(i)-dble(k))/(dble(k)-rd(j)) enddo enddo enddo c do i=0,5 c write(6,*) 'cova',i c do k=1,nk c write(6,*) (cova(k,j,i),j=1,nk) c enddo c enddo return end ********************************************************************** * * * SUBROUTINE ogsy * * * ********************************************************************** subroutine OGSY(cova, runphi, runthe,s,covy) ********************************************************************** c Compute elements of sigma_11 and sigma_22 c Input: cova: array of size (I, I,0:nw) containing covariance of MVARFIMA(0 c runphi: array of size (I, I,0:np) containing current phi values c runthe: array of size (I, I, 0:nq) containing current theta values c Output: s: array of size (I,I,1-np-n:nq-1)) containing eta(k) c covy: array of size(I,I,0:npqmx-1) containing acf of MVARFIMA(p,d, c NOTE: Assumed that runphi(0) and runthe(0) equal -Identity matrix implicit integer(i-n) integer trunc,l,m,o c Note: r equals max(np,nq), nw>= n-1+np+trunc c parameter (npqmx=1,nmx=132) parameter (nkmx=3,nwmx=500) parameter (ntruncmx=350) real*8 runphi(nkmx,nkmx,0:npqmx),runthe(nkmx,nkmx,0:npqmx) real*8 cova(nkmx,nkmx,0:nwmx),covy(nkmx,nkmx,0:npqmx-1), & s(nkmx,nkmx,1-npqmx-nmx:npqmx-1) real*8 A(npqmx*nkmx,npqmx*nkmx),H(npqmx*nkmx,nkmx),det(2), & Ainv(npqmx*nkmx,nkmx),work(nkmx*npqmx) real*8 psi(nkmx,nkmx,0:ntruncmx),w1(nkmx,nkmx),w2(nkmx,nkmx) real*8 w4(nkmx,nkmx) real*8 w3(nkmx,nkmx),c(nkmx,nkmx),d(nkmx,nkmx),zero(nkmx,nkmx), & sum(nkmx,nkmx) integer ipvt(npqmx*nkmx) common/nconst/np,nq,nk,ixd,n,iter1,nrep1,nloop,nrsig,nw call scaler(nk,0.0d0,zero) trunc=ntruncmx c Build A do i = 1, np do j = 1, np+nq if (((j-i) .lt. 0) .or. ((j-i) .gt. nq))then do l=1,nk do m=1,nk A((i-1)*nk+l,(j-1)*nk+m) = 0.0d0 enddo enddo else do l=1,nk do m=1,nk A((i-1)*nk+l,(j-1)*nk+m) =-runthe(l,m,j-i) enddo enddo endif enddo enddo c do i = np+1,np+nq do j = 1, np+nq if (((i-j) .lt. 0) .or. ((i-j) .gt. np))then do l=1,nk do m=1,nk A((i-1)*nk+l,(j-1)*nk+m) = 0.0d0 enddo enddo else do l=1,nk do m=1,nk A((i-1)*nk+l,(j-1)*nk+m) =-runphi(l,m,i-j) enddo enddo endif enddo enddo c c c compute inverse of A call dgefa(A,npqmx*nkmx,(np+nq)*nk,ipvt,info) call dgedi(A,npqmx*nkmx,(np+nq)*nk,ipvt,det,work,01) c c Call mvarmawt (Needed to compute covy) call mvarmawt(trunc,runphi,runthe,psi) c compute covy do k=0,np-1 do ik=1,nk do jk=1,nk sum(ik,jk)=zero(ik,jk) enddo enddo c SUM=ZERO DO I=0,TRUNC DO J=0,TRUNC DO L=1,nk DO M=1,nk W1(L,M)=PSI(L,M,I) W2(L,M)=PSI(M,L,J) IF ((I-J-K) .GE. 0) THEN W3(L,M)=COVA(L,M,I-J-K) ELSE W3(L,M)=COVA(M,L,J+K-I) ENDIF ENDDO ENDDO CALL MATMUL(W1,W3,W4,nk,nk,nk) CALL MATMUL(W4,W2,W1,nk,nk,nk) do ik=1,nk do jk=1,nk sum(ik,jk)=sum(ik,jk)+w1(ik,jk) enddo enddo c SUM=SUM+W1 ENDDO ENDDO do l=1,nk do m=1,nk covy(l,m,k)=sum(l,m) enddo enddo enddo c c compute part of H that contains Ck do k = (1-np)-n, -n c initialize Ck do ik=1,nk do jk=1,nk c(ik,jk)=zero(ik,jk) enddo enddo c c=zero c compute sum of the product runphi*covy do jj = 0, np c Compute necessary covy lag (k+jj) lag=iabs(k+jj) do ik=1,nk do jk=1,nk sum(ik,jk)=zero(ik,jk) enddo enddo c SUM=ZERO DO I=0,TRUNC DO J=0,TRUNC DO L=1,nk DO M=1,nk W1(L,M)=PSI(L,M,I) W2(L,M)=PSI(M,L,J) IF ((I-J-lag) .GE. 0) THEN W3(L,M)=COVA(L,M,I-J-lag) ELSE W3(L,M)=COVA(M,L,J+lag-I) ENDIF ENDDO ENDDO CALL MATMUL(W1,W3,W4,nk,nk,nk) CALL MATMUL(W4,W2,W1,nk,nk,nk) do ik=1,nk do jk=1,nk sum(ik,jk)=sum(ik,jk)+w1(ik,jk) enddo enddo c SUM=SUM+W1 ENDDO ENDDO c compute product runphi(jj)*covy(k+jj) if ((k+jj).lt. 0) then do l=1,nk do m=1,nk do o=1,nk c(l,m)=c(l,m)+runphi(l,o,jj)*sum(m,o) enddo enddo enddo else do l=1,nk do m=1,nk do o=1,nk c(l,m)=c(l,m)+runphi(l,o,jj)*sum(o,m) enddo enddo enddo endif c placing the computed value Ck in H do l=1,nk do m=1,nk H((n+np-1+k)*nk+l,m) = -c(l,m) enddo enddo enddo enddo c c c compute part of H that contains D do k = 1-n, nq-n c initialize Dk do ik=1,nk do jk=1,nk d(ik,jk)=zero(ik,jk) enddo enddo c d = zero c compute sum of the product runthe*cova do i = 0,nq c d = d + runthe(i)*cova'(i-k) if ((k-i) .lt. 0) then do l=1,nk do m=1,nk do o=1,nk d(l,m)=d(l,m)+runthe(l,o,i)*cova(m,o,i-k) enddo enddo enddo else c d = d + runthe(i)*cova(k-i) do l=1,nk do m=1,nk do o=1,nk d(l,m)=d(l,m)+runthe(l,o,i)*cova(o,m,k-i) enddo enddo enddo endif enddo c placing the computed value Dk in H do l=1,nk do m=1,nk H((n+np-1+k)*nk+l, m) = -d(l,m) enddo enddo enddo c compute the first (p+q) elements of s call MATMUL(A,H,Ainv,(np+nq)*nk,(np+nq)*nk,nk) do k=1-np-n,nq-n do l=1,nk do m=1,nk s(l,m,k)=Ainv((k-1+np+n)*nk+l,m) write(6,*) l,m,k,s(l,m,k) enddo enddo enddo c c compute the rest of the elements of s do k = nq-n+1,nq-1 c compute needed Dk c initialize Dk do ik=1,nk do jk=1,nk d(ik,jk)=zero(ik,jk) enddo enddo c d = zero c compute sum of the product runthe*cova do i = 0, nq if ((k-i) .lt. 0) then c d = d + runthe(i)*cova'(i-k) do l=1,nk do m=1,nk do o=1,nk d(l,m)=d(l,m)+runthe(l,o,i)*cova(m,o,i-k) enddo enddo enddo else c d = d + runthe(i)*cova(k-i) do l=1,nk do m=1,nk do o=1,nk d(l,m)=d(l,m)+runthe(l,o,i)*cova(o,m,k-i) enddo enddo enddo endif enddo c compute the sum of the product eta_k and runphi do ik=1,nk do jk=1,nk w3(ik,jk)=zero(ik,jk) enddo enddo c w3=zero do jj = 1, np c s(k) = s(k) + runphi(jj)*s(k-jj)- d do l=1,nk do m=1,nk do o=1,nk w3(l,m)=w3(l,m)+runphi(l,o,jj)*s(o,m,k-jj) enddo enddo enddo enddo do l=1,nk do m=1,nk s(l,m,k)= w3(l,m)-d(l,m) enddo enddo c enddo return end ********************************************************************** * * * SUBROUTINE gsy * * * ********************************************************************** subroutine GSY(cova,runphi,runthe,s,covy) ********************************************************************** c Compute elements of sigma_11 and sigma_22 c Input: cova: array of size (I, I,0:nw) containing covariance of MVARFIMA(0 c runphi: array of size (I, I,0:np) containing current phi values c runthe: array of size (I, I, 0:nq) containing current theta values c Output: s: array of size (I,I,1-np-n:nq-1)) containing eta(k) c covy: array of size(I,I,0:npqmx-1) containing acf of MVARFIMA(p,d, c NOTE: Assumed that runphi(0) and runthe(0) equal -Identity matrix c Note: r equals max(np,nq), nw>= n-1+np+trunc c parameter (npqmx=1,nmx=132) parameter (nkmx=3,nwmx=500) parameter (ntruncmx=350) implicit real*8 (a-h,o-z) integer trunc,l,m,o real*8 runphi(nkmx,nkmx,0:npqmx),runthe(nkmx,nkmx,0:npqmx) real*8 cova(nkmx,nkmx,0:nwmx),covy(nkmx,nkmx,0:npqmx-1), & s(nkmx,nkmx,1-npqmx-nmx:npqmx-1) real*8 psi(nkmx,nkmx,0:ntruncmx),w1(nkmx,nkmx),w2(nkmx,nkmx) real*8 w4(nkmx,nkmx) real*8 w3(nkmx,nkmx),d(nkmx,nkmx),zero(nkmx,nkmx), & sum(nkmx,nkmx) common/nconst/np,nq,nk,ixd,n,iter1,nrep1,nloop,nrsig,nw call scaler(nk,0.0d0,zero) trunc=ntruncmx c c Call mvarmawt (Needed to compute covy) call mvarmawt(trunc,runphi,runthe,psi) c compute covy do k=0,np-1 do ik=1,nk do jk=1,nk sum(ik,jk)=zero(ik,jk) enddo enddo c SUM=ZERO DO I=0,TRUNC DO J=0,TRUNC DO L=1,nk DO M=1,nk W1(L,M)=PSI(L,M,I) W2(L,M)=PSI(M,L,J) IF ((I-J-K) .GE. 0) THEN W3(L,M)=COVA(L,M,I-J-K) ELSE W3(L,M)=COVA(M,L,J+K-I) ENDIF ENDDO ENDDO CALL MATMUL(W1,W3,W4,nk,nk,nk) CALL MATMUL(W4,W2,W1,nk,nk,nk) do ik=1,nk do jk=1,nk sum(ik,jk)=sum(ik,jk)+w1(ik,jk) enddo enddo c SUM=SUM+W1 ENDDO ENDDO c write(6,*) 'in gsy:covy',k do l=1,nk do m=1,nk covy(l,m,k)=sum(l,m) enddo c write(6,*) (covy(l,m,k),m=1,nk) enddo enddo c c approximating method to obtain first (p+q) elements of s s=0.0d0 do k=1-np-n,nq-n do j=0,TRUNC DO L=1,nk DO M=1,nk W1(L,M)=PSI(L,M,J) IF ((k-j) .GE. 0) THEN W2(L,M)=COVA(L,M,k-j) ELSE W2(L,M)=COVA(M,L,J-k) ENDIF ENDDO ENDDO CALL MATMUL(W1,W2,W3,nk,nk,nk) do ik=1,nk do jk=1,nk s(ik,jk,k)=s(ik,jk,k)+w3(ik,jk) enddo enddo enddo c write(6,*) 'in gsy: s ', k c do ik=1,nk c write(6,100) (s(ik,jk,k),jk=1,nk) 100 format(3F10.7) c enddo enddo c c compute the rest of the elements of s do k = nq-n+1,nq-1 c compute needed Dk c initialize Dk do ik=1,nk do jk=1,nk d(ik,jk)=zero(ik,jk) enddo enddo c d = zero c compute sum of the product runthe*cova do i = 0, nq if ((k-i) .lt. 0) then c d = d + runthe(i)*cova'(i-k) do l=1,nk do m=1,nk do o=1,nk d(l,m)=d(l,m)+runthe(l,o,i)*cova(m,o,i-k) enddo enddo enddo else c d = d + runthe(i)*cova(k-i) do l=1,nk do m=1,nk do o=1,nk d(l,m)=d(l,m)+runthe(l,o,i)*cova(o,m,k-i) enddo enddo enddo endif enddo c compute the sum of the product eta_k and runphi do ik=1,nk do jk=1,nk w3(ik,jk)=zero(ik,jk) enddo enddo c w3=zero do jj = 1, np c s(k) = s(k) + runphi(jj)*s(k-jj)- d do l=1,nk do m=1,nk do o=1,nk w3(l,m)=w3(l,m)+runphi(l,o,jj)*s(o,m,k-jj) enddo enddo enddo enddo do l=1,nk do m=1,nk s(l,m,k)= w3(l,m)-d(l,m) enddo enddo c enddo return end ********************************************************************** * * * SUBROUTINE musig * * * ********************************************************************** subroutine musig(runphi,runthe,runy0,runeps0,data,runeps) ********************************************************************** c use:form the error terms: runeps from arfima(0,d,0) process of length c 1-nq:n c parameters of routine: c runphi *input * array of size (nkmx,nkmx,0:npqmx) c runthe *input * array of size (nkmx,nkmx,0:npqmx) c runy0 *input * array of size (nkmx,1-npqmx:0) c runeps0 *input * array of size (nkmx,1-npqmx:0) c data *input * array of size (nkmx,nmx) c runeps *output* array of size (nkmx,1-npqmx:nmx) c parameter (npqmx=1,nmx=132) parameter (nnqmx=132) parameter (nkmx=3) implicit real*8 (a-h,o-z) real*8 data(nkmx,nmx) real*8 runphi(nkmx,nkmx,0:npqmx),runthe(nkmx,nkmx,0:npqmx) real*8 runy0(nkmx,1-npqmx:0),runeps0(nkmx,1-npqmx:0) real*8 runeps(nkmx,1-npqmx:nmx),runy(nkmx,1-npqmx:nmx) real*8 w1(nkmx,nkmx),w2(nkmx),w3(nkmx),s1(nkmx),s3(nkmx) common/nconst/np,nq,nk,ixd,n,iter1,nrep1,nloop,nrsig,nw if (np.ne.0) then do i=1-np,0 do j=1,nk runy(j,i)=runy0(j,i) enddo enddo endif c if (nq.ne.0) then do i=1-nq,0 do j=1,nk runeps(j,i)=runeps0(j,i) enddo enddo endif c do i=1,n do j=1,nk runy(j,i)=data(j,i) enddo enddo c do mt=1,n do j=1,nk s1(j)=runy(j,mt) s3(j)=0.0d0 enddo c if (np.ne.0) then do i=1,np do j=1,nk w2(j)=runy(j,mt-i) do k=1,nk w1(j,k)=runphi(j,k,i) enddo enddo call vmatmul(w1,w2,w3,nk,nk) do ik=1,nk s1(ik)=s1(ik)-w3(ik) enddo enddo endif c if (nq.ne.0) then do i=1,nq do j=1,nk w2(j)=runeps(j,mt-i) do k=1,nk w1(j,k)=runthe(j,k,i) enddo enddo call vmatmul(w1,w2,w3,nk,nk) do ik=1,nk s3(ik)=s3(ik)+w3(ik) enddo enddo endif c do j=1,nk runeps(j,mt)=s1(j)+s3(j) enddo enddo c return end ********************************************************************* * SUBROUTINE grubin ********************************************************************** subroutine grubin(temp1,temp2,gr) * use: compute the Gelman- Rubin statistic * parameter (nrepx=8) implicit real*8 (a-h,o-z) real*8 temp1(nrepx),temp2(nrepx),si(nrepx) common/nconst/np,nq,nk,ixd,n,iter1,nrep1,nloop,nrsig,nw common /grconst/itergr,idf c do jjrep=1,nrep1 c write(6,*)'temp1 in grubin', temp1(jjrep) c write(6,*)'temp2 in grubin', temp2(jjrep) c enddo c form replication means do jjrep=1,nrep1 temp1(jjrep)=temp1(jjrep)/float(itergr) c write(6,*)'temp1, mean',temp1(jjrep) enddo c form overall mean grmean=0.0d0 do jjrep=1,nrep1 grmean=grmean+temp1(jjrep) enddo grmean=grmean/float(nrep1) c write(6,*)'grmean',grmean c form b b=0.0d0 do jjrep=1,nrep1 b=b+(temp1(jjrep)-grmean)**2 enddo b=b*float(itergr)/(float(nrep1)-1.) c write(6,*)'b',b c form replication S_i^2 do jjrep=1,nrep1 si(jjrep)=temp2(jjrep)/float(itergr) - temp1(jjrep)**2 c write(6,*)'si',si(jjrep) enddo c form w w=0.0d0 do jjrep=1,nrep1 w=w+si(jjrep) enddo w=w/float(nrep1) c write(6,*)'w',w c form gr a1=(float(itergr)-1.)/float(itergr) a2=(float(nrep1)+1.)/(float(nrep1*itergr)) a3=b/w a4=float(idf)/(float(idf)-2.) gr=dsqrt((a1+a2*a3)*a4) c write(6,*)'idf in grubin',idf c write(6,*)'gr in grubin',gr return end ********************************************************************** * * * SUBROUTINE perturb * * * ********************************************************************** subroutine perturb (mphi1,mthe1,y01,eps01,sigma1,mean1,rd1, / rsig,rsigmu,omgch,sigfac,msig1l,mrhol,stdsig, / stdrho,imnfl,isigfl) ***************************************************************** * 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,nk,nk,np) * * mthe1 *input/output* array of dim. (iter1,nrep1,nk,nk,nq) * * y01 *input/output* array of dim. (iter1,nrep1,mk,1-np:0) * * eps01 *input/output* array of dim. (iter1,nrep1,nk,1-nq:0) * * sigma1 *input/output* array of dim. (iter1,nrep1,nk,nk) * * mean1 *input/output* array of dim. (iter1,nrep1) * * rsig *input * array of dim. (npq*nk**2,npq*nk**2) * * Cholesky dec. of var(phi,theta) * * omgch *input * array of len. (npq*nk,npq*nk) * * Cholesky dec. of var(eps0,yo) * * other routines used: checksi * ***************************************************************** parameter (npqmx=1,nmx=132) parameter (nkmx=3) parameter (nrsgmx=12) parameter (nrepx=8) implicit real*8 (a-h,o-z) integer imnfl,isigfl real*8 msig1l (1,nrepx,nkmx) real*8 mrhol(1,nrepx) real*8 stdsig(nkmx) real*8 sigma1(1,nrepx,nkmx,nkmx) real*8 newsig(nkmx,nkmx),sigfac(nkmx,nkmx) real*8 mphi1(1,nrepx,nkmx,nkmx,npqmx) real*8 mthe1(1,nrepx,nkmx,nkmx,npqmx) real*8 mean1(1,nrepx,nkmx),rd1(1,nrepx,nkmx) real*8 y01(1,nrepx,nkmx,1-npqmx:0) real*8 eps01(1,nrepx,nkmx,1-npqmx:0) real*8 tempp(nkmx,nkmx,0:npqmx),tempq(nkmx,nkmx,0:npqmx) real*8 rsig(nrsgmx,nrsgmx),tempd(nkmx) real*8 r1(nrepx,nrsgmx) real*8 omgch(npqmx*nkmx,npqmx*nkmx) real*8 r(nrepx,npqmx*nkmx) real*8 rsigmu(nkmx,nkmx) real*8 r2(nrepx,nkmx),stdrho,r4(nkmx),r5,chkident(nkmx,2*nkmx) real*8 sv(2*nkmx),u(nkmx,2*nkmx),v(nkmx,2*nkmx) common/nconst/np,nq,nk,ixd,n,iter1,nrep1,nloop,nrsig,nw c data overdis/1.0d0/ data overdis2/2.0d0/ c indicat1=0 indicat2=0 indicat3=0 indicat4=0 if (np.ne.0) then do j=0,np do ik=1,nk do jk=1,nk tempp(ik,jk,j)=mphi1(1,1,ik,jk,j) enddo enddo enddo endif c if(nq.ne.0) then do j=0,nq do ik=1,nk do jk=1,nk tempq(ik,jk,j)=mthe1(1,1,ik,jk,j) enddo enddo enddo endif if (ixd .eq.1) then do ik=1,nk tempd(ik)=rd1(1,1,ik) enddo endif c c perturb initial values for replications 2,...,nrep1 do 10 irep=2,nrep1 ntry=0 2000 call DRNMVN(1,(np+nq)*nk,omgch, & (npqmx*nkmx),r,nrepx) call DRNMVN(1,nrsig,rsig,nrsgmx,r1,nrepx) if (imnfl .ne.0) & call DRNMVN(1,nk,rsigmu,nkmx,r2,nrepx) ntry=ntry+1 c perturb epso1 if (nq.ne.0) then do j = 1,nq do ik=1,nk eps01(1,irep,ik,j-nq)=eps01(1,1,ik,j-nq)+ 1 r(1,(np*nk)+(j-1)*nk+ik)*overdis enddo enddo endif c perturb y01 if (np.ne.0) then do j = 1,np do ik=1,nk y01(1,irep,ik,j-np)= y01(1,1,ik,j-np)+ 1 r(1,(j-1)*nk+ik)*overdis enddo enddo endif ********************************************************************** c perturb phi, testing for stationarity ********************************************************************** icount=0 if (np.ne.0) then do j=1,np do ik=1,nk do jk=1,nk icount=icount+1 tempp(ik,jk,j)=mphi1(1,1,ik,jk,j)+ 1 r1(irep,icount+jk)*overdis2 enddo enddo enddo c test for stationarity call checkis(tempp,np,indicat1) endif ********************************************************************** c perturb theta, testing for invertibility ********************************************************************** if(nq.ne.0) then do j = 1,nq do ik=1,nk do jk = 1,nk icount=icount+1 tempq(ik,jk,j)=mthe1(1,1,ik,jk,j)+ 1 r1(irep,icount+jk)*overdis2 enddo enddo enddo c test for invertibility call checkis(tempq,nq,indicat2) endif if (ixd.eq.1) then do i=1,nk icount=icount+1 tempd(ik)=rd1(1,1,ik)+r1(irep,icount)*overdis2 if (dabs(tempd(ik)) .ge. .5d0) indicat3=1 enddo endif *********************************************************************** c test for identifiability of ARMA model if ((np+nq) .ne. 0) then if ((np .ne. 0) .and. (nq .eq.0)) then ifl=1 do i=1,nk do j=1,nk chkident(i,j)=tempp(i,j,np) enddo enddo else if ((np .eq. 0) .and. (nq .ne.0)) then ifl=1 do i=1,nk do j=1,nk chkident(i,j)=tempq(i,j,nq) enddo enddo else ifl=2 do i=1,nk do j=1,nk chkident(i,j)=tempp(i,j,np) chkident(i,nk+j)=tempq(i,j,nq) enddo enddo endif endif call DLSVRR(nk,ifl*nk,chkident,nkmx,11,10.0d0*DMACH(4),irank, & sv,u,nkmx,v,nkmx) if (irank .lt. nk) indicat4=1 endif if ((indicat1.eq.1) .or.(indicat2.eq.1) .or. (indicat3.eq.1) & .or. (indicat4.eq.1) ) then if (ntry .lt.50) then overdis2=1.0d0 goto 2000 else write(6,*) 'In perturb: ntry is ', ntry overdis2=2.00d0 ntry=0 goto 2000 endif endif c store perturbed values of phi in mphi1 if (np.ne.0) then do j =1,np do ik=1,nk do jk=1,nk mphi1(1,irep,ik,jk,j)=tempp(ik,jk,j) c write(6,*)'mphi1in pertrub',mphi1(1,irep,ik,jk,j) enddo enddo enddo endif c store perturbed values of theta in mthe1 if (nq.ne.0) then do j =1,nq do ik=1,nk do jk=1,nk mthe1(1,irep,ik,jk,j)=tempq(ik,jk,j) c write(6,*)'mthe1 in perturb',mthe1(1,irep,ik,jk,j) enddo enddo enddo endif c store perturbed values of d in rd1 if(ixd.eq.1) then do ik=1,nk rd1(1,irep,ik)=tempd(ik) enddo endif if (imnfl .ne.0) then c perturb mean1 do ik=1,nk mean1(1,irep,ik)=mean1(1,1,ik)+r2(1,ik)*overdis enddo else do ik=1,nk mean1(1,irep,ik)=0.0d0 enddo endif ************************************************************************ if (isigfl.eq.0) then * perturb sigma1 call drawwish(sigfac,newsig) do ik=1,nk do jk=1,nk sigma1(1,irep,ik,jk)=newsig(ik,jk) enddo enddo else * perturb msig1l do ik=1,nk r4(ik)=drnnof()*stdsig(ik) msig1l(1,irep,ik)=msig1l(1,1,ik)+r4(ik)*overdis enddo * perturb mrhol r5=drnnof()*stdrho mrhol(1,irep)=mrhol(1,1)+r5*overdis endif ************************************************************************ c 10 continue return end ********************************************************************** * * * SUBROUTINE updatemu * * * ********************************************************************** subroutine updatemu (fu,phi,vi,runmn,runphi,runthe,runy0, 1 runeps0,runeps,rsigmu,data,covy,s,a1,a2, 1 det1,de2,desig) ********************************************************************** * 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 * * det1 *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 * * detsig *input * scalar. determinant of wn covariance sigma * * rund *input * array of length nkmx, frac. diff. pars. * 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=1,nmx=132) parameter (nnqmx=132,nwmx=500) parameter (nkmx=3) implicit real*8 (a-h,o-z) real*8 data(nkmx,nmx) real*8 runphi(nkmx,nkmx,0:npqmx),runthe(nkmx,nkmx,0:npqmx) real*8 runy0(nkmx,1-npqmx:0),runeps0(nkmx,1-npqmx:0) real*8 runeps(nkmx,1-npqmx:nmx),runepst(nkmx,1-npqmx:nmx) real*8 runmn(nkmx) real*8 rsigmu(nkmx,nkmx) real*8 vi(nkmx,nkmx,0:nnqmx-1) real*8 phi(nkmx,nkmx,0:nnqmx-1,0:nnqmx-1) real*8 s(nkmx,nkmx,1-npqmx-nmx:npqmx-1) real*8 covy(nkmx,nkmx,0:npqmx-1) real*8 r(1,nkmx) common/nconst/np,nq,nk,ixd,n,iter1,nrep1,nloop,nrsig,nw c do 1000 loop=1,nloop call DRNMVN(1,nk,rsigmu,nkmx,r,1) ********************************************************************** c perturb mu by adding this generation from nMVnormal with 0 mean ********************************************************************** do ik=1,nk runmn(ik)=runmn(ik)+r(1,ik) enddo c do i=1,n do ik=1,nk data(ik,i)=data(ik,i)-r(1,ik) enddo enddo c if (np.ne.0) then do j=1,np do ik=1,nk runy0(ik,j-np)=runy0(ik,j-np)-r(1,ik) enddo enddo endif c call musig(runphi,runthe,runy0,runeps0,data,runepst) call g1(runepst,phi,vi,a1t) if (np.ne.0) call g2(phi,vi,runepst,runy0,covy,s,a2t,de2t) call fcal (fv,a1t,a2t,det1,de2t,desig) c fuv=fv-fu ulg=dlog(drnunf()) c write(6,*)'in updatemu loop,fu,fv,fuv,ulg',loop,fu,fv,fuv,ulg if(ulg.lt.fuv) then fu=fv do j=1-nq,n do ik=1,nk runeps(ik,j)=runepst(ik,j) enddo enddo a1=a1t if(np.ne.0) then a2=a2t de2=de2t endif else do ik=1,nk runmn(ik)=runmn(ik)-r(1,ik) enddo c do in=1,n do ik=1,nk data(ik,in)=data(ik,in)+r(1,ik) enddo enddo c if (np.ne.0) then do j=1,np do ik=1,nk runy0(ik,j-np)=runy0(ik,j-np)+r(1,ik) enddo enddo endif c endif c 1000 continue return end ********************************************************************** * * * SUBROUTINE updtsig * * * ********************************************************************** subroutine updtsig(fu,phi,vi,runphi,runthe, & runy0,runeps,runsig,sigfac,rund,cova,covy,s, & a1,a2,det1,de2,desig) ********************************************************************** * use: generate sample of sigmas by Metropolis algorithm * * parameters of routine: * * fu *input/output* scalar. log of posterior density * * runphi *input * array of length 0:np. * * runthe *input * array of length 0:nq * * rund *input * array of lenght nk ********************************************************************** parameter (npqmx=1,nmx=132) parameter (nkmx=3,nnqmx=132,nwmx=500) implicit real*8 (a-h,o-z) real*8 runsig(nkmx,nkmx),tsig(nkmx,nkmx),xdesig real*8 rund(nkmx),runy0(nkmx,1-npqmx:0) real*8 sigfac(nkmx,nkmx),desig,oldsig(nkmx,nkmx) real*8 runphi(nkmx,nkmx,0:npqmx),runthe(nkmx,nkmx,0:npqmx) real*8 det(2),work(nkmx),runeps(nkmx,1-npqmx:nmx) real*8 covat(nkmx,nkmx,0:nwmx),vit(nkmx,nkmx,0:nnqmx-1) real*8 phit(nkmx,nkmx,0:nnqmx-1,0:nnqmx-1) real*8 st(nkmx,nkmx,1-npqmx-nmx:npqmx-1) real*8 covyt(nkmx,nkmx,0:npqmx-1) real*8 cova (nkmx,nkmx,0:nwmx),vi (nkmx,nkmx,0:nnqmx-1) real*8 phi (nkmx,nkmx,0:nnqmx-1,0:nnqmx-1) real*8 s (nkmx,nkmx,1-npqmx-nmx:npqmx-1) real*8 covy (nkmx,nkmx,0:npqmx-1) integer ipvt(nkmx) common/nconst/np,nq,nk,ixd,n,iter1,nrep1,nloop,nrsig,nw do 1000 loop=1,nloop do ik=1,nk do jk=1,nk oldsig(ik,jk)=runsig(ik,jk) enddo enddo xdesig=desig call drawwish(sigfac,runsig) tsig=runsig call DGEFA(tsig,nkmx,nk,ipvt,info) call DGEDI(tsig,nkmx,nk,ipvt,det,work,10) desig=dlog(det(1)*10.0d0**det(2)) call gcova(covat,rund,runsig) call gc(covat,det1t,vit,phit) call g1(runeps,phit,vit,a1t) c call g1test(covat,runeps,det1t,a1t) if (np.ne.0) then call gsy(covat,runphi,runthe,st,covyt) call g2(phit,vit,runeps,runy0,covyt,st,a2t,de2t) endif call fcal(fv,a1t,a2t,det1t,de2t,desig) fuv=fv-fu ulg=dlog(DRNUNF()) c if(ulg.lt.fuv) then fu=fv a1=a1t det1=det1t c copy phit,vit,covat,covyt,st into corresp. qts. minus t do ik=1,nk do jk=1,nk do i=0,n+nq-1 do j= 0,n+nq-1 phi(ik,jk,i,j)=phit(ik,jk,i,j) enddo enddo enddo enddo c do ik=1,nk do jk=1,nk do i=0,n+nq-1 vi(ik,jk,i)=vit(ik,jk,i) enddo enddo enddo c if (np.ne.0) then a2=a2t de2=de2t if (ixd.eq.1) then do ik=1,nk do jk=1,nk do i=0,nw cova(ik,jk,i)=covat(ik,jk,i) enddo enddo enddo endif c do ik=1,nk do jk=1,nk do i=0,np+nq-1 covy(ik,jk,i)=covyt(ik,jk,i) enddo enddo enddo c do ik=1,nk do jk=1,nk do i=1-np+nq-n,np+nq-1 s(ik,jk,i)= st(ik,jk,i) enddo enddo enddo endif c else runsig=oldsig desig=xdesig endif c 1000 continue return end ********************************************************************** * * * SUBROUTINE updtsig2 * * * ********************************************************************** subroutine updtsig2(fu,phi,vi,runphi,runthe,stdsig, 1 runy0,runeps,runsigl,rund,runrho,cova,covy,s,a1,a2, 1 det1,de2,desig) ********************************************************************** * use: generate sample of sigmas 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 * * det1 *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 * * detsig *input * scalar. determinant of wn covariance sigma * * rund *input * array of length nkmx, frac. diff. pars. * 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=1,nmx=132) parameter (nnqmx=132,nwmx=500) parameter (nkmx=3) implicit real*8 (a-h,o-z) real*8 runsigl(nkmx),runrho real*8 runsig(nkmx,nkmx) real*8 stdsig(nkmx) real*8 runphi(nkmx,nkmx,0:npqmx),runthe(nkmx,nkmx,0:npqmx) real*8 runy0(nkmx,1-npqmx:0) real*8 runeps(nkmx,1-npqmx:nmx) real*8 rund(nkmx) real*8 covat(nkmx,nkmx,0:nwmx),vit(nkmx,nkmx,0:nnqmx-1) real*8 phit(nkmx,nkmx,0:nnqmx-1,0:nnqmx-1) real*8 st(nkmx,nkmx,1-npqmx-nmx:npqmx-1) real*8 covyt(nkmx,nkmx,0:npqmx-1) real*8 cova (nkmx,nkmx,0:nwmx),vi (nkmx,nkmx,0:nnqmx-1) real*8 phi (nkmx,nkmx,0:nnqmx-1,0:nnqmx-1) real*8 s (nkmx,nkmx,1-npqmx-nmx:npqmx-1) real*8 covy (nkmx,nkmx,0:npqmx-1),work(nkmx) real*8 rn(nkmx),tsig(nkmx,nkmx),oldsig(nkmx,nkmx),det(2) integer ipvt(nkmx) common/nconst/np,nq,nk,ixd,n,iter1,nrep1,nloop,nrsig,nw c do 1000 loop=1,nloop xdesig=desig oldsig=runsig do ik=1,nk rn(ik)=drnnof()*stdsig(ik) runsigl(ik)=runsigl(ik)+rn(ik) enddo c rho=(dexp(runrho))/(1.0d0+dexp(runrho)) do ik=1,nk do jk=1,nk t=(dexp(runsigl(ik)))*(dexp(runsigl(jk))) if(ik.eq.jk) runsig(ik,jk)=t if(ik.ne.jk)runsig(ik,jk)=t*rho enddo enddo tsig=runsig call DGEFA(tsig,nkmx,nk,ipvt,info) call DGEDI(tsig,nkmx,nk,ipvt,det,work,10) desig=dlog(det(1)*10.0d0**det(2)) call gcova(covat,rund,runsig) call gc(covat,det1t,vit,phit) call g1 (runeps,phit,vit,a1t) c call g1test(covat,runeps,det1t,a1t) if (np.ne.0) then call gsy(covat,runphi,runthe,st,covyt) call g2(phit,vit,runeps,runy0,covyt,st,a2t,de2t) endif call fcal(fv,a1t,a2t,det1t,de2t,desig) c fuv=fv-fu ulg=dlog(drnunf()) if(ulg.lt.fuv) then fu=fv a1=a1t det1=det1t c copy phit,vit,a2t,de2t,covat,covyt,st into corresp. qts. minus t do ik=1,nk do jk=1,nk do i=0,n+nq-1 do j= 0,n+nq-1 phi(ik,jk,i,j)=phit(ik,jk,i,j) enddo enddo enddo enddo c do ik=1,nk do jk=1,nk do i=0,n+nq-1 vi(ik,jk,i)=vit(ik,jk,i) enddo enddo enddo c if (np.ne.0) then a2=a2t de2=de2t if (ixd.eq.1) then do ik=1,nk do jk=1,nk do i=0,nw cova(ik,jk,i)=covat(ik,jk,i) enddo enddo enddo endif c do ik=1,nk do jk=1,nk do i=0,np+nq-1 covy(ik,jk,i)=covyt(ik,jk,i) enddo enddo enddo c do ik=1,nk do jk=1,nk do i=1-np+nq-n,np+nq-1 s(ik,jk,i)= st(ik,jk,i) enddo enddo enddo endif else do ik=1,nk runsigl(ik)=runsigl(ik)-rn(ik) enddo runsig=oldsig desig=xdesig c endif c 1000 continue return end ********************************************************************** * * * SUBROUTINE updtrho * * * ********************************************************************** subroutine updtrho(fu,phi,vi,runphi,runthe,stdrho, 1 runy0,runeps,runsigl,runrho,rund,cova,covy,s,a1,a2, 1 det1,de2,desig) ********************************************************************** * use: generate sample of sigmas 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 * * det1 *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 * * detsig *input * scalar. determinant of wn covariance sigma * * rund *input * array of length nkmx, frac. diff. pars. * 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=1,nmx=132) parameter (nnqmx=132,nwmx=500) parameter (nkmx=3) implicit real*8 (a-h,o-z) real*8 runsigl(nkmx),runrho,work(nkmx) real*8 runsig(nkmx,nkmx),oldsig(nkmx,nkmx),tsig(nkmx,nkmx),det(2) real*8 runphi(nkmx,nkmx,0:npqmx),runthe(nkmx,nkmx,0:npqmx) real*8 runy0(nkmx,1-npqmx:0) real*8 runeps(nkmx,1-npqmx:nmx) real*8 rund(nkmx) real*8 rn real*8 cova(nkmx,nkmx,0:nwmx),vi(nkmx,nkmx,0:nnqmx-1) real*8 phi(nkmx,nkmx,0:nnqmx-1,0:nnqmx-1) real*8 s(nkmx,nkmx,1-npqmx-nmx:npqmx-1) real*8 covy(nkmx,nkmx,0:npqmx-1) real*8 covat(nkmx,nkmx,0:nwmx),vit(nkmx,nkmx,0:nnqmx-1) real*8 phit(nkmx,nkmx,0:nnqmx-1,0:nnqmx-1) real*8 st(nkmx,nkmx,1-npqmx-nmx:npqmx-1) real*8 covyt(nkmx,nkmx,0:npqmx-1) integer ipvt(nkmx) common/nconst/np,nq,nk,ixd,n,iter1,nrep1,nloop,nrsig,nw do 1000 loop=1,nloop xdesig=desig oldsig=runsig rn=drnnof()*stdrho runrho = runrho +rn c rho=(dexp(runrho))/(1.0d0+dexp(runrho)) do ik=1,nk do jk=1,nk t=(dexp(runsigl(ik)))*(dexp(runsigl(jk))) if(ik.eq.jk) runsig(ik,jk)=t if(ik.ne.jk)runsig(ik,jk)=t*rho enddo enddo tsig=runsig call DGEFA(tsig,nkmx,nk,ipvt,info) call DGEDI(tsig,nkmx,nk,ipvt,det,work,10) desig=dlog(det(1)*10.0d0**det(2)) call gcova(covat,rund,runsig) call gc(covat,det1t,vit,phit) call g1 (runeps,phit,vit,a1t) c call g1test(covat,runeps,det1t,a1t) if (np.ne.0) then call gsy(covat,runphi,runthe,st,covyt) call g2(phit,vit,runeps,runy0,covyt,st,a2t,de2t) endif call fcal(fv,a1t,a2t,det1t,de2t,desig) c fuv=fv-fu ulg=dlog(drnunf()) c write(6,*)'in updtrho loop,fu,fv,fuv,ulg',loop,fu,fv,fuv,ulg if(ulg.lt.fuv) then fu=fv a1=a1t det1=det1t c copy phit,vit,a2t,de2t,covat,covyt,st into corresp. qts. minus t do ik=1,nk do jk=1,nk do i=0,n+nq-1 do j= 0,n+nq-1 phi(ik,jk,i,j)=phit(ik,jk,i,j) enddo enddo enddo enddo c do ik=1,nk do jk=1,nk do i=0,n+nq-1 vi(ik,jk,i)=vit(ik,jk,i) enddo enddo enddo c if (np.ne.0) then a2=a2t de2=de2t if (ixd.eq.1) then do ik=1,nk do jk=1,nk do i=0,nw cova(ik,jk,i)=covat(ik,jk,i) enddo enddo enddo endif c do ik=1,nk do jk=1,nk do i=0,np+nq-1 covy(ik,jk,i)=covyt(ik,jk,i) enddo enddo enddo c do ik=1,nk do jk=1,nk do i=1-np+nq-n,np+nq-1 s(ik,jk,i)= st(ik,jk,i) enddo enddo enddo endif else runrho=runrho-rn runsig=oldsig desig=xdesig endif c 1000 continue return end ********************************************************************** * * * SUBROUTINE update * * * ********************************************************************** subroutine update(fu,phi,vi,runphi,runthe,runy0,runeps0, / runeps,runsig,rsig,data,rund,cova,covy,s, / a1,a2,det1,de2,desig) ********************************************************************** * use: generate sample for (PHI,THETA,D) vectors * * usig metropolis algorithm * * parameters of routine: * * fu *input/output* scalar. log of posterior density * * phi *input * array of dimension (I,I,0:n+nq-1,0:n+nq-1) contains * forward partial correlation matrices needed for * sigma 22 inverse * v *input * array of dimension (I,I,0:n+nq-1) contains block * diagonal elements needed for sigma 22 inverse * det1 *input/output * scalar. contains determinant term in g1 * * function * * rd *input/output* vector. generated value of fractional * * difference parameter vector 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 * * detsig *input * scalar. determinant of WN covariance sigma * * 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=1,nmx=132) parameter (nnqmx=132,nwmx=500) parameter (nkmx=3) parameter (nrsgmx=12) implicit real*8 (a-h,o-z) real*8 data(nkmx,nmx) real*8 runsig(nkmx,nkmx) real*8 runphi(nkmx,nkmx,0:npqmx),runthe(nkmx,nkmx,0:npqmx) real*8 runy0(nkmx,1-npqmx:0),runeps0(nkmx,1-npqmx:0) real*8 runeps(nkmx,1-npqmx:nmx),runepst(nkmx,1-npqmx:nmx) real*8 rund(nkmx) real*8 cova(nkmx,nkmx,0:nwmx),vi(nkmx,nkmx,0:nnqmx-1) real*8 phi(nkmx,nkmx,0:nnqmx-1,0:nnqmx-1) real*8 s(nkmx,nkmx,1-npqmx-nmx:npqmx-1) real*8 covy(nkmx,nkmx,0:npqmx-1) real*8 r(1,nrsgmx) real*8 tempp(nkmx,nkmx,npqmx) real*8 tempq(nkmx,nkmx,npqmx) real*8 tempd(nkmx) real*8 rsig(nrsgmx,nrsgmx) real*8 covat(nkmx,nkmx,0:nwmx), covyt(nkmx,nkmx,0:npqmx-1) real*8 st(nkmx,nkmx ,1-npqmx-nmx:npqmx-1) real*8 phit(nkmx,nkmx,0:nnqmx-1,0:nnqmx-1) real*8 vit(nkmx,nkmx,0:nnqmx-1), chkident(nkmx,2*nkmx) real*8 sv(2*nkmx),u(nkmx,2*nkmx),v(nkmx,2*nkmx) common/nconst/np,nq,nk,ixd,n,iter1,nrep1,nloop,nrsig,nw ********************************************************************** c ntry=0 loop=1 2000 call DRNMVN(1,nrsig,rsig,nrsgmx,r,1) indicat1=0 indicat2=0 indicat3=0 indicat4=0 ntry=ntry+1 if(np.ne.0) then do i = 1,np do ik=1,nk do jk=1,nk tempp(ik,jk,i)=runphi(ik,jk,i) enddo enddo enddo endif c if(nq.ne.0) then do i = 1,nq do ik=1,nk do jk=1,nk tempq(ik,jk,i)=runthe(ik,jk,i) enddo enddo enddo endif c if (ixd.eq.1) then do ik=1,nk tempd(ik)=rund(ik) enddo endif c ********************************************************************** c perturb phi, testing for stationarity ********************************************************************** icount=0 if (np.ne.0) then do j=1,np do ik=1,nk do jk=1,nk icount=icount+1 runphi(ik,jk,j)=runphi(ik,jk,j)+r(1,icount) enddo enddo enddo call checkis(runphi,np,indicat1) endif ********************************************************************** c perturb theta, testing for invertibility ********************************************************************** 5650 if (nq.ne.0) then do j=1,nq do ik=1,nk do jk=1,nk icount=icount+1 runthe(ik,jk,j)=runthe(ik,jk,j)+r(1,icount) enddo enddo enddo call checkis(runthe,nq,indicat2) endif *********************************************************************** c perturb rund1 *********************************************************************i if (ixd.eq.1) then do ik=1,nk rund (ik)=rund(ik)+r(1,icount+ik) enddo c write(6,*) 'rund: ', (rund(ik),ik=1,nk) *********************************************************************** c test for stationarity of each scalar rund ********************************************************************** do ik=1,nk if(dabs(rund(ik)).ge.0.5d0) indicat3=1 enddo endif *********************************************************************** c test for identifiability of ARMA model if ((np+nq) .ne. 0) then if ((np .ne. 0) .and. (nq .eq.0)) then ifl=1 do i=1,nk do j=1,nk chkident(i,j)=runphi(i,j,np) enddo enddo else if ((np .eq. 0) .and. (nq .ne.0)) then ifl=1 do i=1,nk do j=1,nk chkident(i,j)=runthe(i,j,nq) enddo enddo else ifl=2 do i=1,nk do j=1,nk chkident(i,j)=runphi(i,j,np) chkident(i,nk+j)=runthe(i,j,nq) enddo enddo endif endif call DLSVRR(nk,ifl*nk,chkident,nkmx,11,10.0d0*DMACH(4),irank, & sv,u,nkmx,v,nkmx) if (irank .lt. nk) indicat4=1 endif c c write(6,*) 'i1,i2,i3,i4: ', indicat1,indicat2,indicat3,indicat4 if ((indicat1.eq.1) .or.(indicat2.eq.1) .or. (indicat3.eq.1) .or. & (indicat4.eq.1)) then do j=1,np do ik=1,nk do jk=1,nk runphi(ik,jk,j)=tempp(ik,jk,j) enddo enddo enddo do j=1,nq do ik=1,nk do jk=1,nk runthe(ik,jk,j)=tempq(ik,jk,j) enddo enddo enddo do ik=1,nk rund(ik)=tempd(ik) enddo if (ntry .lt. 50) then goto 2000 else go to 5656 endif endif 5656 if (ntry.ge.50) then write(6,*)'ntry ',ntry goto 3000 endif 3000 call musig(runphi,runthe,runy0,runeps0,data,runepst) call gcova(covat,rund,runsig) call gc(covat,det1t,vit,phit) call g1(runepst,phit,vit,a1t) c call g1test(covat,runepst,det1t,a1t) if (np.ne.0) then call gsy(covat,runphi,runthe,st,covyt) call g2(phit,vit,runepst,runy0,covyt,st,a2t,de2t) endif call fcal (fv,a1t,a2t,det1t,de2t,desig) c fuv=fv-fu ulg=dlog(drnunf()) c write(6,*)'in update:loop,fu,fv,a1t,det1t',loop,fu,fv,a1t,det1t if(ulg.lt.fuv) then fu=fv do j=1-np-nq,n do ik=1,nk runeps(ik,j)=runepst(ik,j) enddo enddo a1=a1t det1=det1t c copy phit,vit,a2t,de2t,covat,covyt,st into corresp. qts. minus t do ik=1,nk do jk=1,nk do i=0,n+nq-1 do j= 0,n+nq-1 phi(ik,jk,i,j)=phit(ik,jk,i,j) enddo enddo enddo enddo c do ik=1,nk do jk=1,nk do i=0,n+nq-1 vi(ik,jk,i)=vit(ik,jk,i) enddo enddo enddo c if (np.ne.0) then a2=a2t de2=de2t if (ixd.eq.1) then do ik=1,nk do jk=1,nk do i=0,nw cova(ik,jk,i)=covat(ik,jk,i) enddo enddo enddo endif c do ik=1,nk do jk=1,nk do i=0,np+nq-1 covy(ik,jk,i)=covyt(ik,jk,i) enddo enddo enddo c do ik=1,nk do jk=1,nk do i=1-np-nq-n,np+nq-1 s(ik,jk,i)= st(ik,jk,i) enddo enddo enddo endif else if (np.ne.0) then do j=1,np do ik=1,nk do jk=1,nk runphi(ik,jk,j)=tempp(ik,jk,j) enddo enddo enddo endif c if (nq.ne.0) then do j=1,nq do ik=1,nk do jk=1,nk runthe(ik,jk,j)=tempq(ik,jk,j) enddo enddo enddo endif c if (ixd.eq.1) then do ik=1,nk rund(ik)=tempd(ik) enddo endif c endif loop=loop+1 if((loop.lt.nloop) .and.(ntry .lt. 50)) goto 2000 c return end ********************************************************************** * * * SUBROUTINE upeps * * * ********************************************************************** subroutine upeps(fu,phi,vi,runphi,runthe,runy0,runeps0, 1 omgch,data,runeps,covy,s,a1,a2, 1 det1,de2,desig) ********************************************************************** * use: generate sample for latent variables y0 and eps0 using * * metropolis algorithm * * parameters of routine: * * fu *input/output* scalar. log of posterior density * * vi *input-array of size (I,I,0:n+q-1) containing v_j^{-1} * * phi *input- array of size (I,I,0:n+nq-1,0:n+nq-1) containing * * *forward partial correlation matrices * * *(needed for sigma_{2,2}^-1 * * det1 *input * scalar. contains ln(det) term in g1 function* * runphi *input * array of dimension (nkmx,nkmx,0:npqmx) * * runthe *input * array of dimension (nkmx,nkmx,0:npqmx) * * runy0 *input/output* array of dimension (nkmx,1-npqmx:0) * * * generated sample for y0 vector of latent vars * * runeps0 *input/output* array of dimension (nkmx,1-npqmx:0) * * * generated sample for eps0 vector of latent vars * * omgch *input * array of dimension (npqkmx*(npqkmx+1)/2. * * * contains estimated Cholesky decomposition of covar. * * for sampling (y0,eps0) * * data *input * array of dimension (nkmx,nmx) * * runeps *input/output* array of dimension (nkmx,1-npqmx:nmx) * * detsig *input * scalar. determinant of covariance SIGMA * * rund *input * array of dimension nk. frac. diff. pars. * * cova *input * array of dimension (nkmx,nkmx,0:nw) contains * * * acf from mvarfima(0,d,0) process * * covy *input * array of dimension nkmx,nkmx,0:nnqmx-1) contains* * acf of mvarfima(p,d,q) process * * s *input * array of dim. (nkmx,nkmx,1-npqmx-nmx:npqmx-1) * * * other routines used: rand,dmvnor,musig,g1,g2,fcal * ********************************************************************** parameter (npqmx=1,nmx=132) parameter (nnqmx=132) parameter (nwmx=500) parameter (nkmx=3) implicit real*8 (a-h,o-z) real*8 data(nkmx,nmx) real*8 runphi(nkmx,nkmx,0:npqmx),runthe(nkmx,nkmx,0:npqmx) real*8 runy0(nkmx,1-npqmx:0),runeps0(nkmx,1-npqmx:0) real*8 runeps(nkmx,1-npqmx:nmx),runepst(nkmx,1-npqmx:nmx) real*8 omgch(npqmx*nkmx,nkmx*npqmx) real*8 vi(nkmx,nkmx,0:nnqmx-1) real*8 phi(nkmx,nkmx,0:nnqmx-1,0:nnqmx-1) real*8 s(nkmx,nkmx,1-npqmx-nmx:npqmx-1) real*8 covy(nkmx,nkmx,0:npqmx-1) real*8 r(1,npqmx*nkmx) common/nconst/np,nq,nk,ixd,n,iter1,nrep1,nloop,nrsig,nw npk=np*nk do 50 loop =1,nloop call DRNMVN(1,(np+nq)*nk,omgch, & nkmx*npqmx,r,1) if(nq.ne.0) then do j=1,nq do ik=1,nk runeps0(ik,j-nq)=runeps0(ik,j-nq)+ 1 r(1,(np*nk)+(j-1)*nk+ik) enddo enddo endif c if(np.ne.0) then do j = 1,np do ik=1,nk runy0(ik,j-np)=runy0(ik,j-np)+ 1 r(1,(j-1)*nk+ik) enddo enddo endif c call musig(runphi,runthe,runy0,runeps0,data,runepst) call g1 (runepst,phi,vi,a1t) if(np.ne.0) call g2(phi,vi,runepst,runy0,covy,s,a2t,de2t) call fcal (fv,a1t,a2t,det1,de2t,desig) fuv=fv-fu ulg=dlog(drnunf()) c write(6,*)'in upeps loop,fu,fv,fuv,ulg',loop,fu,fv,fuv,ulg c if (ulg.lt.fuv)then fu=fv do j=1-nq,n do ik=1,nk runeps(ik,j)=runepst(ik,j) enddo enddo a1=a1t if(np.ne.0) then a2=a2t de2=de2t endif else if(nq.ne.0) then do j=1,nq do ik=1,nk runeps0(ik,j-nq)=runeps0(ik,j-nq)- 1 r(1,np*(nk**2)+(j-1)*nk+ik) enddo enddo endif c if(np.ne.0) then do j=1,np do ik=1,nk runy0(ik,j-np)=runy0(ik,j-np)- 1 r(1,(j-1)*nk+ik) enddo enddo endif endif c 50 continue c return end ********************************************************************** * * * SUBROUTINE xomega * * * ********************************************************************** subroutine XOMEGA(omega, runphi,runthe,rd,sigma) c Compute var-cov of latent variables eps0 and y0 c Input: runphi: array of size (I, I,0:np) containing initial phi values c runthe: array of size (I, I, 0:nq) containing initial theta values c rd : vector of length I containing initial d values c sigma: var-cov matrix of errors c Output: omega: vector of length ((npq*I)*((npq+1)*I)/2) containing lower c triangular elements of omega in order (1,1),(2,1),(2,2), etc c Other routines used: MATMUL (matrix multiplication) c DGEDI (matrix inverse/uses DAXPY,DSCAL,DSWAP) c DGEFA (matrix determinant/used DAXPY,DSCAL,IDAMAX) c MVARMAWT(MA wts of MVARMA process) c SCALER c NOTE: Assumed that runphi(0) and runthe(0) equal -Identity matrix implicit integer (i-n) parameter (npqmx=1,nmx=132) parameter (nnqmx=132) parameter (nwmx=500) parameter (nkmx=3) parameter (ntruncmx=350) integer trunc real*8 runphi(nkmx,nkmx,0:npqmx),runthe(nkmx,nkmx,0:npqmx) real*8 cova(nkmx,nkmx,0:nwmx) real*8 omega(npqmx*nkmx, npqmx*nkmx) real*8 covy(nkmx,nkmx,0:npqmx-1) real*8 s(nkmx,nkmx,1-npqmx-nmx:npqmx-1) real*8 sigma(nkmx,nkmx),w4(nkmx,nkmx),sum(nkmx,nkmx) real*8 psi(nkmx,nkmx,0:ntruncmx),rd(nkmx),w1(nkmx,nkmx) real*8 w2(nkmx,nkmx) real*8 w3(nkmx,nkmx),zero(nkmx,nkmx) common/nconst/np,nq,nk,ixd,n,iter1,nrep1,nloop,nrsig,nw trunc=ntruncmx call scaler(nk,0.0d0,zero) c Compute cov matrix for MVARFIMA(0,d_i,0) call gcova(cova,rd,sigma) c Call mvarmawt (Needed to compute covy) call mvarmawt(trunc,runphi,runthe,psi) c compute covy do k=0,np-1 do ik=1,nk do jk=1,nk sum(ik,jk)=zero(ik,jk) enddo enddo c SUM=ZERO DO I=0,TRUNC DO J=0,TRUNC DO L=1,nk DO M=1,nk W1(L,M)=PSI(L,M,I) W2(L,M)=PSI(M,L,J) IF ((I-J-K) .GE. 0) THEN W3(L,M)=COVA(L,M,I-J-K) ELSE W3(L,M)=COVA(M,L,J+K-I) ENDIF ENDDO ENDDO CALL MATMUL(W1,W3,W4,nk,nk,nk) CALL MATMUL(W4,W2,W1,nk,nk,nk) do ik=1,nk do jk=1,nk sum(ik,jk)=sum(ik,jk)+w1(ik,jk) enddo enddo c SUM=SUM+W1 ENDDO ENDDO do l=1,nk do m=1,nk covy(l,m,k)=sum(l,m) enddo enddo enddo c Approximating method to obtain first (p+q) elements of s s=0.0d0 do k=1-np,nq do j=0,TRUNC DO L=1,nk DO M=1,nk W1(L,M)=PSI(L,M,J) IF ((k-j) .GE. 0) THEN W2(L,M)=COVA(L,M,k-j) ELSE W2(L,M)=COVA(M,L,J-k) ENDIF ENDDO ENDDO CALL MATMUL(W1,W2,W3,nk,nk,nk) do ik=1,nk do jk=1,nk s(ik,jk,k)=s(ik,jk,k)+w3(ik,jk) enddo enddo enddo enddo c Place elements in OMEGA do i=1,np do j=1,i do l=1,nk do m=1,nk omega((i-1)*nk+l,(j-1)*nk+m)=covy(l,m,i-j) omega((j-1)*nk+l,(i-1)*nk+m)=covy(m,l,i-j) enddo enddo enddo enddo c sigma_12 (pI*qI) and sigma_21 (qI*pI) do i=1,np do j=1,nq do l=1,nk do m=1,nk omega((i-1)*nk+l,(np+j-1)*nk+m)=s(l,m,i-j-(np-nq)) omega((np+j-1)*nk+l,(i-1)*nk+m)=s(m,l,i-j-(np-nq)) enddo enddo enddo enddo c sigma_22 (qI*qI) do i=1,nq do j=1,i do l=1,nk do m=1,nk omega((np+i-1)*nk+l,(np+j-1)*nk+m)=cova(l,m,i-j) omega((np+j-1)*nk+l,(np+i-1)*nk+m)=cova(m,l,i-j) enddo enddo enddo enddo return end ********************************************************************** SUBROUTINE MVARMAWT(LMAX,RUNPHI,RUNTHE,PSI) ********************************************************************** C parameter (npqmx=1,nmx=132) parameter (nkmx=3) parameter (ntruncmx=350) implicit real*8 (a-h,o-z) integer R, S, I,L,M,O,MXPQ1,LMAX real*8 psi(nkmx,nkmx,0:ntruncmx),w(nkmx,nkmx),zero(nkmx,nkmx) real*8 runphi(nkmx,nkmx,0:npqmx),runthe(nkmx,nkmx,0:npqmx) real*8 phi (nkmx,nkmx,0:npqmx),theta(nkmx,nkmx,0:npqmx) common/nconst/np,nq,nk,ixd,n,iter1,nrep1,nloop,nrsig,nw R=MAX (NP,NQ) c copys runphi and runthe to phi and theta (setting phi(0) c and theta(0)=-I) call scaler(nk,0.0d0,zero) do l=1,nk do m=1,nk do o=0,np phi(l,m,o)=runphi(l,m,o) enddo do o=0,nq theta(l,m,o)=runthe(l,m,o) enddo enddo enddo C C ZEROISE PLANES OF PHI,THETA BEYOND THE INPUT PLANES (IF ANY). C IF (NP .EQ.NQ) GOTO 30 C C NONE IN THIS CASE. C C IF (R .EQ.NQ) CALL LINZER(PHI(1, 1,NP + 1),NK *NK * (R -NP)) if (r .eq.nq) then do i=1,(r-np) do l=1,nk do m=1,nk phi(l,m,np+i)=0.0d0 enddo enddo enddo endif C IF (R .EQ.NP) CALL LINZER(THETA(1, 1,NQ + 1),NK *NK * (R -NQ)) if (r .eq.np) then do i=1,(r-nq) do l=1,nk do m=1,nk theta(l,m,nq+i)=0.0d0 enddo enddo enddo endif C C NEXT, FORM THE PSI-MATRICES, UP TO THE LMAX'TH. C 30 CALL SCALER(NK,1.0d0,w) do l=1,nk do m=1,nk psi(l,m,0)=w(l,m) enddo enddo C PSI(0) = IDENTITY C MXPQ1=MAX(np,nq+1) DO S = 1, LMAX C C FORM PSI(S) C do ik=1,nk do jk=1,nk w(ik,jk)=zero(ik,jk) enddo enddo c W=ZERO DO I=1,MIN(S,NP) DO L=1,NK DO M=1,NK DO O=1,NK W(L,M)=W(L,M)+PHI(L,O,I)*PSI(O,M,S-I) ENDDO ENDDO ENDDO IF (S .LT. MXPQ1) THEN DO L=1,NK DO M=1,NK PSI(L,M,S)=W(L,M)-THETA(L,M,S) ENDDO ENDDO ELSE DO L=1,NK DO M=1,NK PSI(L,M,S)=W(L,M) ENDDO ENDDO ENDIF ENDDO ENDDO C 130 RETURN END C *********************************************************************** subroutine prx0(grphi1,grphi2,grthe1,grthe2,grd1,grd2,gry01,gry02, & greps01,greps02, grmn1,grmn2, grsig1,grsig2,grsigl1,grsigl2, & grrho1,grrho2,imnfl,isigfl) *********************************************************************** 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 mean1 *input* array of dimension (iter1,nrep1) c sigma1 *input* array of dimension (iter1,nrep1,nk,nk) c c____ print the samplers ____________________________________________* parameter (npqmx=1,nmx=132) parameter (nkmx=3) parameter (nrepx=8) implicit real*8 (a-h,o-z) real*8 sumphi1(nkmx,nkmx,npqmx),sumphi2(nkmx,nkmx,npqmx) real*8 sumthe1(nkmx,nkmx,npqmx),sumthe2(nkmx,nkmx,npqmx) real*8 summn1(nkmx),summn2(nkmx) real*8 sumsig1(nkmx,nkmx),sumsigl1(nkmx) real*8 sumsig2(nkmx,nkmx), sumsigl2(nkmx) real*8 sumrho1,sumrho2,sumy01(nkmx,npqmx),sumy02(nkmx,npqmx) real*8 sumeps01(nkmx,npqmx),sumeps02(nkmx,npqmx) real*8 sumd1(nkmx),sumd2(nkmx) real*8 xmean(nkmx,nkmx),xvar(nkmx,nkmx) real*8 rmean(nkmx),rvar