! Fotran 90 code for the estimating the SCSAR model used in the empirical study of ! "A Social Interactions Model with Endogenous Friendship Formation and Selectivity" ! forthcoming in Journal of Applied Econometrics. ! ! This f90 file should be complied along with the other two module files -- ! random.f90 contributed by Alan Miller, and tool.f90 ! ! ! Written by ! Chih-Sheng Hsieh ! Department of Economics ! The Chinese University of Hong Kong ! 11/17/2014 ! ! Any question can be directed to cshsieh@cuhk.edu.hk ! !***************************************************************************************** include 'mkl_vsl.f90' PROGRAM SCSAR USE random USE TOOL USE mkl95_LAPACK USE MKL_VSL_TYPE USE MKL_VSL USE omp_lib IMPLICIT NONE !****************** GENERAL VARIABLES ****************! CHARACTER(150) :: path1='D:/research/social_interaction_model_with_selectivity/' CHARACTER(150) :: path2='/empirical/data/GPA_50_with_missing/' CHARACTER(150) :: path3='/empirical/result/missing/D1/' CHARACTER(40) :: numchr INTEGER, PARAMETER :: GG=73, TT=150000, MAX=141, DD=1, THIN=100, Nthreads=200 REAL, PARAMETER :: pi=3.14159265 REAL :: c_1, c_2, acc_1, acc_2, rand, c_3(GG), acc_3(GG), acc_4, acc_3v INTEGER :: NN, N(GG), k, count, missing_num(GG), missing INTEGER :: brng, method, me, i, j, g, t, d, flag, ier, r, seedsize, seed(2), errorflag, INFO, seed_mkl, t_round INTEGER :: clock_max, clock_rate, clock_start, clock_stop, ithread, iblock,ibegin,iend REAL, DIMENSION(MAX,MAX,GG):: age, sex, race, sport, network REAL, DIMENSION(MAX,19,GG):: group REAL, ALLOCATABLE:: C1(:,:), C2(:,:), C3(:,:), C4(:,:), W(:,:), IPIV(:) REAL, ALLOCATABLE:: Y(:), overage(:) REAL, ALLOCATABLE:: X(:,:) LOGICAL :: first=.true. TYPE (VSL_STREAM_STATE) :: stream INTEGER :: errcode REAL(kind=8) :: det1, det2, pp !*** hyperparameters in the prior distribution ***! REAL, DIMENSION(6+DD) :: gamma_0 REAL, DIMENSION(28) :: beta_0 REAL, DIMENSION(6+DD,6+DD) :: G_0 REAL, DIMENSION(28,28) :: B_0, INV_B_0 REAL, DIMENSION(1+DD) :: sigma_0 REAL, DIMENSION(1+DD,1+DD) :: Sig_0 REAL :: eta_0, A_0, alpha_0 INTEGER :: rho_0 !*** posterior draws of parameters *******! REAL, DIMENSION(6+DD,TT):: gamma_T REAL, DIMENSION(TT) :: lambda_T, acc_rate1, acc_rate2, acc_rate4 REAL, DIMENSION(28, TT) :: beta_T REAL, DIMENSION(DD+1,TT) :: sigma_T REAL, DIMENSION(GG,TT) :: alpha_T, acc_rate3 REAL(KIND=8) :: logp_T(TT) !****** VARIABLE IN M-H FOR Y_T *******************************! REAL :: h_i, Y_CT(MAX,GG), Y_T(MAX,GG) REAL, ALLOCATABLE :: INV_S(:,:), INV_V(:,:) REAL, ALLOCATABLE :: T_i(:), C_i(:), Y_C(:) REAL, ALLOCATABLE :: MU_Y(:), Y_m(:) INTEGER, ALLOCATABLE :: P(:) !****** VARIABLE IN M-H FOR Z_G *******************************! REAL, DIMENSION(MAX,DD,GG):: zz REAL, DIMENSION(MAX,DD,GG,TT/THIN):: zz_T REAL, ALLOCATABLE :: zz_1(:,:), zz_2(:,:), V(:,:), S1(:,:), S2(:,:), XX(:,:) REAL, ALLOCATABLE :: V1(:,:), V2(:,:), Inv_V2(:,:), V3(:,:), Vc(:,:), Inv_Vc(:,:) REAL, ALLOCATABLE :: psi(:,:), psi_1(:,:), psi_2(:,:), zero(:) REAL, ALLOCATABLE :: ep_1(:), ep_2(:) REAL, ALLOCATABLE :: ep_1v(:), ep_2v(:) REAL(kind=8) :: like_Y1, like_Y2, q_1, q_2 !********VARIABLE IN M-H FOR GAMMA & LAMBDA ***************! REAL :: lambda_1, lambda_2 REAL :: gamma_1(6+DD), gamma_2(6+DD), COV_gamma_1(6+DD,6+DD) REAL(kind=8) :: pp_l(GG), pp_G(GG), pp_sig(GG) !*************VARIABLE IN M-H FOR BETA **************! REAL :: XVX(28,28,GG), XVY(28,GG), INV_B(28,28), INV_B2(28,28), beta_temp(28) REAL, ALLOCATABLE :: XX2(:,:), YY(:), S(:,:) !*************VARIABLE IN M-H FOR SIGMA_E & SIGMA_EZ & ALPHA_G**************! REAL :: sigma_1(1+DD), sigma_2(1+DD), COV_sigma_1(1+DD,1+DD) REAL :: R_g REAL, ALLOCATABLE :: VV1(:,:), VV2(:,:) !***********************************************************! seedsize=2 seed(1)=2.1103e+005 seed(2)=2.1103e+005 CALL RANDOM_SEED CALL RANDOM_SEED(SIZE=seedsize) CALL RANDOM_SEED(PUT=seed) seed_mkl=12345 brng=VSL_BRNG_MCG31 method=VSL_METHOD_SGAUSSIANMV_BOXMULLER me=VSL_MATRIX_STORAGE_FULL errcode=vslnewstream(stream,brng,seed_mkl) CALL networksize(GG, MAX, path1, path2, N) CALL define_variable(GG, age, sex, race, sport, network, group, N, path1, path2, MAX) gamma_0=0.0 beta_0=0.0 G_0=eye(6+DD)*10.0 Sig_0=eye(1+DD)*10.0 B_0=eye(28)*10.0 CALL FINDInv(B_0,INV_B_0,28,ier) eta_0=0.1 rho_0=2 zz=0.0 A_0=400.0 alpha_0=0.0 c_3=0.5 acc_3=0.0 sigma_0(1)=1.0 sigma_0(2:DD+1)=0.0 ! initial value of the draws ! gamma_T(:,1)=(/-1.0626 , 0.2787 , 0.6454 , 0.6575 , -0.0101 , 0.0614 , 5.2349/) !DD=1! !gamma_T(:,1)=(/ 0.5981 , 0.2547 , 0.6930 , 0.7413 , -0.0326 , 0.0530 , 3.4313 , 3.1545/) !DD=2! !gamma_T(:,1)=(/ 2.0678 , 0.2630 , 0.7486 , 0.9049 , -0.0462 , 0.0569 , 3.0883 , 2.9994 , 2.8828/) !DD=3! !gamma_T(:,1)=(/ 3.1420 , 0.2535 , 0.7465 , 0.9589 , -0.0462 , 0.0569 , 2.7371 , 2.4012 , 2.1742 , 1.6778/) !DD=4! !gamma_T(:,1)=(/ 3.9384 , 0.2481 , 0.7706 , 1.0170 , 2.9352 , 2.2109 , 1.9746 , 1.4077 , 1.2004/) !DD=5! lambda_T(1)=0.02 sigma_T(1,1)=1.0 sigma_T(2:DD+1,1)=0.0 OPEN (1000, file = trim(adjustl(path1)) // trim(adjustl(path3)) // 'gamma.txt', status='unknown') OPEN (2000, file = trim(adjustl(path1)) // trim(adjustl(path3)) // 'lambda.txt', status='unknown') OPEN (3000, file = trim(adjustl(path1)) // trim(adjustl(path3)) // 'beta.txt', status='unknown') OPEN (4000, file = trim(adjustl(path1)) // trim(adjustl(path3)) // 'sigma.txt', status='unknown') OPEN (5000, file = trim(adjustl(path1)) // trim(adjustl(path3)) // 'latent_Z.txt', status='unknown') OPEN (6000, file = trim(adjustl(path1)) // trim(adjustl(path3)) // 'alpha.txt', status='unknown') OPEN (7000, file = trim(adjustl(path1)) // trim(adjustl(path3)) // 'logp.txt', status='unknown') missing_num=0 DO g=1,GG NN=N(g) ALLOCATE(Y(NN)); count=0 Y=group(1:NN,18,g) DO i=1,NN IF (Y(i)==0) count=count+1 END DO missing_num(g)=count DEALLOCATE(Y); END DO DO t=2, TT call system_clock ( clock_start, clock_rate, clock_max ) DO lambda_1=random_normal()*0.1+lambda_T(t-1) IF (t>2) THEN lambda_2=random_normal()*sqrt(var(lambda_T(1:t-1),t-1))*2.38+lambda_T(t-1) lambda_1=0.95*lambda_2+0.05*lambda_1 END IF IF ((lambda_1>=0.00) .and. (lambda_1<=0.10)) EXIT END DO DO COV_gamma_1=EYE(6+DD)*0.01/(6+DD) call spotrf('U',6+DD, COV_gamma_1, 6+DD, info) errcode=vsrnggaussianmv(method,stream,1,gamma_1, 6+DD, me, gamma_0, COV_gamma_1) gamma_1=gamma_1+gamma_T(:,t-1) IF (t>2*(DD+6)) THEN COV_gamma_1=cov(gamma_T(:,1:t-1),6+DD,t-1)*2.38**2/(6+DD) call spotrf('U',6+DD, COV_gamma_1, 6+DD, info) errcode=vsrnggaussianmv(method,stream,1,gamma_2, 6+DD, me, gamma_0, COV_gamma_1) gamma_2=gamma_2+gamma_T(:,t-1) gamma_1=0.95*gamma_2+0.05*gamma_1 END IF IF (DD==1) THEN IF (gamma_1(6+DD)>=0) EXIT ELSEIF (DD==2) THEN IF (gamma_1(5+DD)>=gamma_1(6+DD)) EXIT ELSEIF (DD==3) THEN IF ((gamma_1(4+DD)>=gamma_1(5+DD)) .and. (gamma_1(5+DD)>=gamma_1(6+DD)) ) EXIT ELSEIF (DD==4) THEN IF ((gamma_1(3+DD)>=gamma_1(4+DD)) .and. (gamma_1(4+DD)>=gamma_1(5+DD)) .and. (gamma_1(5+DD)>=gamma_1(6+DD))) EXIT ELSEIF (DD==5) THEN IF ((gamma_1(2+DD)>=gamma_1(3+DD)) .and. (gamma_1(3+DD)>=gamma_1(4+DD)) .and. (gamma_1(4+DD)>=gamma_1(5+DD)) .and. (gamma_1(5+DD)>=gamma_1(6+DD)) ) EXIT END IF END DO gamma_2=gamma_T(:,t-1) ! THE M-H ALGORITHM FOR SAMPLING Y_missing ! Call OMP_SET_NUM_THREADS(Nthreads) !$omp parallel default(shared) private(g,NN, W, Y, X, XX, zz_1, S1, INV_S, V, INV_V, & !$omp & T_i, C_i, P, MU_Y, Y_C, IPIV, ier, INFO, i, h_i, rand, missing, Vc, Inv_Vc, V1, V2, V3, Inv_V2, Y_m, & !$omp & ithread, iblock, ibegin, iend) ithread = OMP_GET_THREAD_NUM() iblock = (GG+Nthreads-1)/Nthreads ibegin = ithread*iblock+1 iend = min((ithread+1)*iblock,GG) DO g = ibegin,iend NN=N(g) ALLOCATE(W(NN,NN)); ALLOCATE(Y(NN)) ALLOCATE(X(NN,14)); ALLOCATE(XX(NN,28)) ALLOCATE(S1(NN,NN)); ALLOCATE(INV_S(NN,NN)) ALLOCATE(V(NN,NN)); ALLOCATE(INV_V(NN,NN)) ALLOCATE(P(NN)); ALLOCATE(MU_Y(NN)) ALLOCATE(zz_1(NN,DD)) ALLOCATE(IPIV(NN)) IF (missing_num(g)>=1) THEN missing=missing_num(g) ALLOCATE(Vc(missing,missing)); ALLOCATE(Inv_Vc(missing,missing)); ALLOCATE(V1(missing,NN-missing)); ALLOCATE(V2(NN-missing,NN-missing)); ALLOCATE(Inv_V2(NN-missing,NN-missing)) ALLOCATE(V3(NN-missing,missing)); ALLOCATE(Y_m(missing)); ALLOCATE(T_i(missing-1)); ALLOCATE(Y_C(missing-1)); ALLOCATE(C_i(missing-1)) W=network(1:NN,1:NN,g) Y=group(1:NN,18,g) X=group(1:NN,1:14,g) zz_1=zz(1:NN,1:DD,g) S1=EYE(NN)-lambda_T(t-1)*W CALL FINDInv(S1,INV_S,NN,ier) V=matmul(matmul(INV_S,(sigma_T(1,t-1)-sum(sigma_T(2:DD+1,t-1)*sigma_T(2:DD+1,t-1)))*EYE(NN)),transpose(INV_S)) XX(:,1:14)=X XX(:,15:28)=matmul(W,X) MU_Y=matmul(XX,beta_T(:,t-1))+matmul(zz_1,sigma_T(2:DD+1,t-1))+alpha_T(g,t-1) CALL SGESV(NN,1,S1,NN,IPIV,MU_Y,NN,INFO) Vc=V(1:missing,1:missing) V1=V(1:missing,missing+1:NN) V2=V(missing+1:NN,missing+1:NN) V3=V(missing+1:NN,1:missing) CALL FINDInv(V2,Inv_V2,NN-missing,ier) Vc=Vc-matmul(matmul(V1,Inv_V2),V3) Y_m=MU_Y(1:missing)+matmul(matmul(V1,Inv_V2),Y(missing+1:NN)-MU_Y(missing+1:NN)) IF (missing==1) THEN Y_CT(1,g)=sum(sqrt(Vc))*random_normal() ELSE CALL FINDInv(Vc,Inv_Vc,missing,ier) DO i=1, missing IF (i==1) THEN T_i=Inv_Vc(i,i+1:missing) Y_C=Y_CT(i+1:missing,g) END IF IF (i==missing) THEN T_i=Inv_Vc(i,1:missing-1) Y_C=Y_CT(1:missing-1,g) END IF IF((i>1) .and. (i=1) Y(1:missing_num(g))=Y_T(1:missing_num(g),g) acc_3v=0.0 DO k=1, NN DO d=1, DD zz_1(k,d)=random_normal()*c_3(g)+zz_2(k,d) END DO pp=0.0 DO i=1,NN DO j=1,NN IF (j/=i) THEN IF (DD==1) THEN psi_1(i,j)=gamma_T(1,t-1)+gamma_T(2,t-1)*c1(i,j)+gamma_T(3,t-1)*c2(i,j)+gamma_T(4,t-1)*c3(i,j)+gamma_T(5,t-1)*overage(i)+gamma_T(6,t-1)*overage(j) & & -gamma_T(7,t-1)*abs(zz_1(i,1)-zz_1(j,1)) psi_2(i,j)=gamma_T(1,t-1)+gamma_T(2,t-1)*c1(i,j)+gamma_T(3,t-1)*c2(i,j)+gamma_T(4,t-1)*c3(i,j)+gamma_T(5,t-1)*overage(i)+gamma_T(6,t-1)*overage(j) & & -gamma_T(7,t-1)*abs(zz_2(i,1) -zz_2(j,1)) ELSEIF (DD==2) THEN psi_1(i,j)=gamma_T(1,t-1)+gamma_T(2,t-1)*c1(i,j)+gamma_T(3,t-1)*c2(i,j)+gamma_T(4,t-1)*c3(i,j)+gamma_T(5,t-1)*overage(i)+gamma_T(6,t-1)*overage(j) & & -gamma_T(7,t-1)*abs(zz_1(i,1)-zz_1(j,1))-gamma_T(8,t-1)*abs(zz_1(i,2)-zz_1(j,2)) psi_2(i,j)=gamma_T(1,t-1)+gamma_T(2,t-1)*c1(i,j)+gamma_T(3,t-1)*c2(i,j)+gamma_T(4,t-1)*c3(i,j)+gamma_T(5,t-1)*overage(i)+gamma_T(6,t-1)*overage(j) & & -gamma_T(7,t-1)*abs(zz_2(i,1)-zz_2(j,1))-gamma_T(8,t-1)*abs(zz_2(i,2)-zz_2(j,2)) ELSEIF (DD==3) THEN psi_1(i,j)=gamma_T(1,t-1)+gamma_T(2,t-1)*c1(i,j)+gamma_T(3,t-1)*c2(i,j)+gamma_T(4,t-1)*c3(i,j)+gamma_T(5,t-1)*overage(i)+gamma_T(6,t-1)*overage(j) & & -gamma_T(7,t-1)*abs(zz_1(i,1)-zz_1(j,1))-gamma_T(8,t-1)*abs(zz_1(i,2)-zz_1(j,2)) -gamma_T(9,t-1)*abs(zz_1(i,3)-zz_1(j,3)) psi_2(i,j)=gamma_T(1,t-1)+gamma_T(2,t-1)*c1(i,j)+gamma_T(3,t-1)*c2(i,j)+gamma_T(4,t-1)*c3(i,j)+gamma_T(5,t-1)*overage(i)+gamma_T(6,t-1)*overage(j) & & -gamma_T(7,t-1)*abs(zz_2(i,1)-zz_2(j,1))-gamma_T(8,t-1)*abs(zz_2(i,2)-zz_2(j,2)) -gamma_T(9,t-1)*abs(zz_2(i,3)-zz_2(j,3)) ELSEIF (DD==4) THEN psi_1(i,j)=gamma_T(1,t-1)+gamma_T(2,t-1)*c1(i,j)+gamma_T(3,t-1)*c2(i,j)+gamma_T(4,t-1)*c3(i,j)+gamma_T(5,t-1)*overage(i)+gamma_T(6,t-1)*overage(j) & & -gamma_T(7,t-1)*abs(zz_1(i,1)-zz_1(j,1))-gamma_T(8,t-1)*abs(zz_1(i,2)-zz_1(j,2)) -gamma_T(9,t-1)*abs(zz_1(i,3)-zz_1(j,3)) & & -gamma_T(10,t-1)*abs(zz_1(i,4)-zz_1(j,4)) psi_2(i,j)=gamma_T(1,t-1)+gamma_T(2,t-1)*c1(i,j)+gamma_T(3,t-1)*c2(i,j)+gamma_T(4,t-1)*c3(i,j)+gamma_T(5,t-1)*overage(i)+gamma_T(6,t-1)*overage(j) & & -gamma_T(7,t-1)*abs(zz_2(i,1)-zz_2(j,1))-gamma_T(8,t-1)*abs(zz_2(i,2)-zz_2(j,2)) -gamma_T(9,t-1)*abs(zz_2(i,3)-zz_2(j,3)) & & -gamma_T(10,t-1)*abs(zz_2(i,4)-zz_2(j,4)) ELSEIF (DD==5) THEN psi_1(i,j)=gamma_T(1,t-1)+gamma_T(2,t-1)*c1(i,j)+gamma_T(3,t-1)*c2(i,j)+gamma_T(4,t-1)*c3(i,j)+gamma_T(5,t-1)*overage(i)+gamma_T(6,t-1)*overage(j) & & -gamma_T(7,t-1)*abs(zz_1(i,1)-zz_1(j,1))-gamma_T(8,t-1)*abs(zz_1(i,2)-zz_1(j,2)) -gamma_T(9,t-1)*abs(zz_1(i,3)-zz_1(j,3)) & & -gamma_T(10,t-1)*abs(zz_1(i,4)-zz_1(j,4)) -gamma_T(11,t-1)*abs(zz_1(i,5)-zz_1(j,5)) psi_2(i,j)=gamma_T(1,t-1)+gamma_T(2,t-1)*c1(i,j)+gamma_T(3,t-1)*c2(i,j)+gamma_T(4,t-1)*c3(i,j)+gamma_T(5,t-1)*overage(i)+gamma_T(6,t-1)*overage(j) & & -gamma_T(7,t-1)*abs(zz_2(i,1)-zz_2(j,1))-gamma_T(8,t-1)*abs(zz_2(i,2)-zz_2(j,2)) -gamma_T(9,t-1)*abs(zz_2(i,3)-zz_2(j,3)) & & -gamma_T(10,t-1)*abs(zz_2(i,4)-zz_2(j,4)) -gamma_T(11,t-1)*abs(zz_2(i,5)-zz_2(j,5)) ENDIF q_1=psi_1(i,j)*W(i,j)-log(1+exp(psi_1(i,j))) q_2=psi_2(i,j)*W(i,j)-log(1+exp(psi_2(i,j))) pp=pp+(q_1-q_2) ENDIF END DO END DO S1=eye(NN)-lambda_T(t-1)*W XX(:,1:14)=X XX(:,15:28)=matmul(W,X) ep_1=matmul(S1,Y)-matmul(XX,beta_T(:,t-1))-matmul(zz_1,sigma_T(2:DD+1,t-1))-alpha_T(g,t-1) ep_2=matmul(S1,Y)-matmul(XX,beta_T(:,t-1))-matmul(zz_2,sigma_T(2:DD+1,t-1))-alpha_T(g,t-1) ep_1v=ep_1; ep_2v=ep_2 V1=V; V2=V CALL SGESV(NN,1,V1,NN,IPIV,ep_1v,NN,INFO) CALL SGESV(NN,1,V2,NN,IPIV,ep_2v,NN,INFO) like_Y1=-dot_product(ep_1, ep_1v)/2 like_Y2=-dot_product(ep_2, ep_2v)/2 pp=pp+(like_Y1-like_Y2)+logmvnpdf(zz_1(k,:), eye(DD), DD)-logmvnpdf(zz_2(k,:),eye(DD), DD) CALL random_number(rand) IF (log(rand)<=pp) THEN zz_2=zz_1 acc_3v=acc_3v+1 END IF zz_1=zz_2 END DO if (2*acc_3v>NN) acc_3(g)=acc_3(g)+1 acc_rate3(g,t)=acc_3(g)/t if (acc_rate3(g,t)<0.2) c_3(g)=c_3(g)/1.01 if (acc_rate3(g,t)>0.3) c_3(g)=c_3(g)*1.01 zz(1:NN,1:DD,g)=zz_1 t_round=t/THIN IF (t-t_round*THIN==0) THEN zz_T(:,:,g,t/THIN)=zz(:,:,g) END IF !***** THE M-H ALGORITHM FOR SAMPLING GAMMA & LAMBDA ******! pp=0.0 DO i=1,NN DO j=1, NN IF (j/=i) THEN IF (DD==1) THEN psi_1(i,j)=gamma_1(1)+gamma_1(2)*c1(i,j)+gamma_1(3)*c2(i,j)+gamma_1(4)*c3(i,j)+gamma_1(5)*overage(i)+gamma_1(6)*overage(j) & & -gamma_1(7)*abs(zz_1(i,1)-zz_1(j,1)) psi_2(i,j)=gamma_2(1)+gamma_2(2)*c1(i,j)+gamma_2(3)*c2(i,j)+gamma_2(4)*c3(i,j)+gamma_2(5)*overage(i)+gamma_2(6)*overage(j) & & -gamma_2(7)*abs(zz_1(i,1)-zz_1(j,1)) ELSEIF (DD==2) THEN psi_1(i,j)=gamma_1(1)+gamma_1(2)*c1(i,j)+gamma_1(3)*c2(i,j)+gamma_1(4)*c3(i,j)+gamma_1(5)*overage(i)+gamma_1(6)*overage(j) & & -gamma_1(7)*abs(zz_1(i,1)-zz_1(j,1))-gamma_1(8)*abs(zz_1(i,2)-zz_1(j,2)) psi_2(i,j)=gamma_2(1)+gamma_2(2)*c1(i,j)+gamma_2(3)*c2(i,j)+gamma_2(4)*c3(i,j)+gamma_2(5)*overage(i)+gamma_2(6)*overage(j) & & -gamma_2(7)*abs(zz_1(i,1)-zz_1(j,1))-gamma_2(8)*abs(zz_1(i,2)-zz_1(j,2)) ELSEIF (DD==3) THEN psi_1(i,j)=gamma_1(1)+gamma_1(2)*c1(i,j)+gamma_1(3)*c2(i,j)+gamma_1(4)*c3(i,j)+gamma_1(5)*overage(i)+gamma_1(6)*overage(j) & & -gamma_1(7)*abs(zz_1(i,1)-zz_1(j,1))-gamma_1(8)*abs(zz_1(i,2)-zz_1(j,2))-gamma_1(9)*abs(zz_1(i,3)-zz_1(j,3)) psi_2(i,j)=gamma_2(1)+gamma_2(2)*c1(i,j)+gamma_2(3)*c2(i,j)+gamma_2(4)*c3(i,j)+gamma_2(5)*overage(i)+gamma_2(6)*overage(j) & & -gamma_2(7)*abs(zz_1(i,1)-zz_1(j,1))-gamma_2(8)*abs(zz_1(i,2)-zz_1(j,2))-gamma_2(9)*abs(zz_1(i,3)-zz_1(j,3)) ELSEIF (DD==4) THEN psi_1(i,j)=gamma_1(1)+gamma_1(2)*c1(i,j)+gamma_1(3)*c2(i,j)+gamma_1(4)*c3(i,j)+gamma_1(5)*overage(i)+gamma_1(6)*overage(j) & & -gamma_1(7)*abs(zz_1(i,1)-zz_1(j,1))-gamma_1(8)*abs(zz_1(i,2)-zz_1(j,2))-gamma_1(9)*abs(zz_1(i,3)-zz_1(j,3)) & & -gamma_1(10)*abs(zz_1(i,4)-zz_1(j,4)) psi_2(i,j)=gamma_2(1)+gamma_2(2)*c1(i,j)+gamma_2(3)*c2(i,j)+gamma_2(4)*c3(i,j)+gamma_2(5)*overage(i)+gamma_2(6)*overage(j) & & -gamma_2(7)*abs(zz_1(i,1)-zz_1(j,1))-gamma_2(8)*abs(zz_1(i,2)-zz_1(j,2))-gamma_2(9)*abs(zz_1(i,3)-zz_1(j,3)) & & -gamma_2(10)*abs(zz_1(i,4)-zz_1(j,4)) ELSEIF (DD==5) THEN psi_1(i,j)=gamma_1(1)+gamma_1(2)*c1(i,j)+gamma_1(3)*c2(i,j)+gamma_1(4)*c3(i,j)+gamma_1(5)*overage(i)+gamma_1(6)*overage(j) & & -gamma_1(7)*abs(zz_1(i,1)-zz_1(j,1))-gamma_1(8)*abs(zz_1(i,2)-zz_1(j,2))-gamma_1(9)*abs(zz_1(i,3)-zz_1(j,3)) & & -gamma_1(10)*abs(zz_1(i,4)-zz_1(j,4)) -gamma_1(11)*abs(zz_1(i,5)-zz_1(j,5)) psi_2(i,j)=gamma_2(1)+gamma_2(2)*c1(i,j)+gamma_2(3)*c2(i,j)+gamma_2(4)*c3(i,j)+gamma_2(5)*overage(i)+gamma_2(6)*overage(j) & & -gamma_2(7)*abs(zz_1(i,1)-zz_1(j,1))-gamma_2(8)*abs(zz_1(i,2)-zz_1(j,2))-gamma_2(9)*abs(zz_1(i,3)-zz_1(j,3)) & & -gamma_2(10)*abs(zz_1(i,4)-zz_1(j,4))-gamma_2(11)*abs(zz_1(i,5)-zz_1(j,5)) ENDIF q_1=psi_1(i,j)*W(i,j)-log(1+exp(psi_1(i,j))) q_2=psi_2(i,j)*W(i,j)-log(1+exp(psi_2(i,j))) pp=pp+(q_1-q_2) ENDIF END DO END DO pp_G(g)=pp XX(:,1:14)=X XX(:,15:28)=matmul(W,X) S1=eye(NN)-lambda_1*W S2=eye(NN)-lambda_T(t-1)*W ep_1=matmul(S1,Y)-matmul(XX,beta_T(:,t-1))-matmul(zz_1,sigma_T(2:DD+1,t-1))-alpha_T(g,t-1) ep_2=matmul(S2,Y)-matmul(XX,beta_T(:,t-1))-matmul(zz_1,sigma_T(2:DD+1,t-1))-alpha_T(g,t-1) ep_1v=ep_1 ep_2v=ep_2 V1=V V2=V CALL SGESV(NN,1,V1,NN,IPIV,ep_1v,NN,INFO) CALL SGESV(NN,1,V2,NN,IPIV,ep_2v,NN,INFO) like_Y1=log(FindDet(S1, NN))-dot_product(ep_1, ep_1v)/2 like_Y2=log(FindDet(S2, NN))-dot_product(ep_2, ep_2v)/2 pp_l(g)=(like_Y1-like_Y2) DEALLOCATE(V); DEALLOCATE(C1); DEALLOCATE(C2) DEALLOCATE(C3); DEALLOCATE(C4) DEALLOCATE(W); DEALLOCATE(Y) DEALLOCATE(S1); DEALLOCATE(S2) DEALLOCATE(X); DEALLOCATE(XX) DEALLOCATE(zz_1); DEALLOCATE(zz_2) DEALLOCATE(V1); DEALLOCATE(V2) DEALLOCATE(ep_1); DEALLOCATE(ep_2) DEALLOCATE(ep_1v); DEALLOCATE(ep_2v) DEALLOCATE(IPIV); DEALLOCATE(overage) DEALLOCATE(psi_1); DEALLOCATE(psi_2) END DO !$omp end parallel pp=sum(pp_G)+logmvnpdf(gamma_1, G_0, 6+DD)-logmvnpdf(gamma_2, G_0, 6+DD) CALL random_number(rand) IF (log(rand)<=pp) THEN gamma_T(:,t)=gamma_1 acc_1=acc_1+1.0 ELSE gamma_T(:,t)=gamma_T(:,t-1) ENDIF acc_rate1(t)=acc_1/t CALL random_number(rand) IF (log(rand)<=sum(pp_l)) THEN lambda_T(t)=lambda_1 acc_2=acc_2+1.0 ELSE lambda_T(t)=lambda_T(t-1) END IF acc_rate2(t)=acc_2/t !***************GIBBS STEP TO SIMULTE BETA *********************! XVX=0.0 XVY=0.0 Call OMP_SET_NUM_THREADS(Nthreads) !$omp parallel default(shared) private(g,NN, V, Y, YY, & !$omp & X, XX, XX2, W, S, V1, V2, zz_1, IPIV, INFO, & !$omp & ithread, iblock, ibegin, iend) ithread = OMP_GET_THREAD_NUM() iblock = (GG+Nthreads-1)/Nthreads ibegin = ithread*iblock+1 iend = min((ithread+1)*iblock,GG) DO g = ibegin,iend NN=N(g) ALLOCATE(V(NN,NN)) ALLOCATE(W(NN,NN)); ALLOCATE(Y(NN)) ALLOCATE(X(NN,14)); ALLOCATE(YY(NN)) ALLOCATE(S(NN,NN)); ALLOCATE(XX(NN,28)) ALLOCATE(XX2(NN,28)); ALLOCATE(zz_1(NN,DD)) ALLOCATE(V1(NN,NN)); ALLOCATE(V2(NN,NN)) ALLOCATE(IPIV(NN)) V=(sigma_T(1,t-1)-sum(sigma_T(2:DD+1,t-1)*sigma_T(2:DD+1,t-1)))*EYE(NN) W=network(1:NN,1:NN,g) Y=group(1:NN,18,g) X=group(1:NN,1:14,g) zz_1=zz(1:NN,1:DD,g) IF (missing_num(g)>=1) Y(1:missing_num(g))=Y_T(1:missing_num(g),g) S=eye(NN)-lambda_T(t)*W XX(:,1:14)=X XX(:,15:28)=matmul(W,X) XX2=XX V1=V CALL SGESV(NN,28,V1,NN,IPIV,XX2,NN,INFO) YY=matmul(S,Y)-matmul(zz_1,sigma_T(2:DD+1,t-1))-alpha_T(g,t-1) V2=V CALL SGESV(NN,1,V2,NN,IPIV,YY,NN,INFO) XVX(:,:,g)=matmul(transpose(XX) , XX2) XVY(:,g)=matmul(transpose(XX) , YY) DEALLOCATE(V) DEALLOCATE(W); DEALLOCATE(Y) DEALLOCATE(X); DEALLOCATE(YY) DEALLOCATE(S); DEALLOCATE(XX) DEALLOCATE(XX2); DEALLOCATE(zz_1) DEALLOCATE(V1); DEALLOCATE(V2) DEALLOCATE(IPIV) END DO !$omp end parallel CALL FINDInv(sum(XVX,3)+INV_B_0,INV_B,28,ier) INV_B2=INV_B CALL spotrf('U',28, INV_B2, 28, info) errcode=vsrnggaussianmv(method, stream, 1, beta_temp, 28, me, beta_0, INV_B2) beta_T(:,t)=beta_temp+matmul(INV_B,sum(XVY,2)) !******************M-H STEP FOR SIGMA_E & SIGMA_EZ *****************! ALLOCATE(zero(1+DD)) zero=0.0 DO COV_sigma_1=EYE(1+DD)*0.01/(1+DD) call spotrf('U',1+DD, COV_sigma_1, 1+DD, info) errcode=vsrnggaussianmv(method,stream,1,sigma_1, 1+DD, me, zero, COV_sigma_1) sigma_1=sigma_1+sigma_T(:,t-1) IF (t>2*(1+DD)) THEN COV_sigma_1=cov(sigma_T(:,1:t-1),1+DD,t-1)*2.38**2/(1+DD) call spotrf('U',1+DD, COV_sigma_1, 1+DD, info) errcode=vsrnggaussianmv(method,stream, 1, sigma_2, 1+DD, me, zero, COV_sigma_1) sigma_2=sigma_2+sigma_T(:,t-1) sigma_1=0.95*sigma_2+0.05*sigma_1 END IF IF (DD==1) THEN IF( (sigma_1(1)>=sum(sigma_1(2:DD+1)*sigma_1(2:DD+1))) .and. (sigma_1(2)>=0) ) EXIT ELSEIF (DD==2) THEN IF( (sigma_1(1)>=sum(sigma_1(2:DD+1)*sigma_1(2:DD+1))) .and. (sigma_1(2)>=0) .and. (sigma_1(3)>=0) ) EXIT ELSEIF (DD==3) THEN IF( (sigma_1(1)>=sum(sigma_1(2:DD+1)*sigma_1(2:DD+1))) .and. (sigma_1(2)>=0) .and. (sigma_1(3)>=0) .and. (sigma_1(4)>=0) ) EXIT ELSEIF (DD==4) THEN IF( (sigma_1(1)>=sum(sigma_1(2:DD+1)*sigma_1(2:DD+1))) .and. (sigma_1(2)>=0) .and. (sigma_1(3)>=0) .and. (sigma_1(4)>=0) .and. (sigma_1(5)>=0) ) EXIT ELSEIF (DD==5) THEN IF( (sigma_1(1)>=sum(sigma_1(2:DD+1)*sigma_1(2:DD+1))) .and. (sigma_1(2)>=0) .and. (sigma_1(3)>=0) .and. (sigma_1(4)>=0) .and. (sigma_1(5)>=0) .and. (sigma_1(6)>=0) ) EXIT END IF END DO DEALLOCATE(zero) Call OMP_SET_NUM_THREADS(Nthreads) !$omp parallel default(shared) private(g,NN, S, V1, V2, Y, VV1, VV2, & !$omp & zz_1, X, XX, W, ep_1, ep_2, ep_1v, ep_2v, IPIV, INFO, like_Y1, like_Y2, & !$omp & ithread, iblock, ibegin, iend) ithread = OMP_GET_THREAD_NUM() iblock = (GG+Nthreads-1)/Nthreads ibegin = ithread*iblock+1 iend = min((ithread+1)*iblock,GG) DO g = ibegin,iend NN=N(g) ALLOCATE(V1(NN,NN)); ALLOCATE(V2(NN,NN)) ALLOCATE(VV1(NN,NN)); ALLOCATE(VV2(NN,NN)) ALLOCATE(W(NN,NN)); ALLOCATE(Y(NN)) ALLOCATE(X(NN,14)); ALLOCATE(XX(NN,28)) ALLOCATE(S(NN,NN)); ALLOCATE(zz_1(NN,DD)) ALLOCATE(ep_1(NN)); ALLOCATE(ep_2(NN)) ALLOCATE(ep_1v(NN)); ALLOCATE(ep_2v(NN)) ALLOCATE(IPIV(NN)) V1=(sigma_1(1)-sum(sigma_1(2:DD+1)*sigma_1(2:DD+1)))*EYE(NN) V2=(sigma_T(1,t-1)-sum(sigma_T(2:DD+1,t-1)*sigma_T(2:DD+1,t-1)))*EYE(NN) W=network(1:NN,1:NN,g) Y=group(1:NN,18,g) X=group(1:NN,1:14,g) zz_1=zz(1:NN,1:DD,g) XX(:,1:14)=X XX(:,15:28)=matmul(W,X) IF (missing_num(g)>=1) Y(1:missing_num(g))=Y_T(1:missing_num(g),g) S=eye(NN)-lambda_T(t)*W ep_1=matmul(S,Y)-matmul(XX,beta_T(:,t))-matmul(zz_1,sigma_1(2:DD+1))-alpha_T(g,t-1) ep_2=matmul(S,Y)-matmul(XX,beta_T(:,t))-matmul(zz_1,sigma_T(2:DD+1,t-1))-alpha_T(g,t-1) ep_1v=ep_1 ep_2v=ep_2 VV1=V1 VV2=V2 CALL SGESV(NN,1,VV1,NN,IPIV,ep_1v,NN,INFO) CALL SGESV(NN,1,VV2,NN,IPIV,ep_2v,NN,INFO) like_Y1=(-0.5)*NN*log(sigma_1(1)-sum(sigma_1(2:DD+1)*sigma_1(2:DD+1)))-dot_product(ep_1, ep_1v)/2 like_Y2=(-0.5)*NN*log(sigma_T(1,t-1)-sum(sigma_T(2:DD+1,t-1)*sigma_T(2:DD+1,t-1)))-dot_product(ep_2, ep_2v)/2 pp_sig(g)=(like_Y1-like_Y2) DEALLOCATE(V1); DEALLOCATE(V2) DEALLOCATE(VV1); DEALLOCATE(VV2) DEALLOCATE(W); DEALLOCATE(Y) DEALLOCATE(X); DEALLOCATE(XX) DEALLOCATE(S); DEALLOCATE(zz_1) DEALLOCATE(ep_1); DEALLOCATE(ep_2) DEALLOCATE(ep_1v);DEALLOCATE(ep_2v) DEALLOCATE(IPIV) END DO !$omp end parallel pp=sum(pp_sig)+logmvnpdf( sigma_1-sigma_0, Sig_0, 1+DD)-logmvnpdf(sigma_T(1:1+DD,t-1)-sigma_0, Sig_0, 1+DD) CALL random_number(rand) IF (log(rand)<=pp) THEN sigma_T(:,t)=sigma_1 acc_4=acc_4+1.0 ELSE sigma_T(:,t)=sigma_T(:,t-1) ENDIF acc_rate4(t)=acc_4/t !*************GIBBS STEP FOR ALPHA_G*********************! Call OMP_SET_NUM_THREADS(Nthreads) !$omp parallel default(shared) private(g,NN, Y, YY, & !$omp & X, XX, W, S, zz_1, IPIV, R_g, & !$omp & ithread, iblock, ibegin, iend) ithread = OMP_GET_THREAD_NUM() iblock = (GG+Nthreads-1)/Nthreads ibegin = ithread*iblock+1 iend = min((ithread+1)*iblock,GG) DO g = ibegin,iend NN=N(G) ALLOCATE(W(NN,NN)); ALLOCATE(Y(NN)) ALLOCATE(X(NN,14)); ALLOCATE(XX(NN,28)) ALLOCATE(S(NN,NN)); ALLOCATE(YY(NN)) ALLOCATE(zz_1(NN,DD)) R_g=(A_0**(-1)+(sigma_T(1,t)-sum(sigma_T(2:DD+1,t)*sigma_T(2:DD+1,t)))**(-1.0)*NN)**(-1.0) W=network(1:NN,1:NN,g) Y=group(1:NN,18,g) X=group(1:NN,1:14,g) zz_1=zz(1:NN,1:DD,g) IF (missing_num(g)>=1) Y(1:missing_num(g))=Y_T(1:missing_num(g),g) S=eye(NN)-lambda_T(t)*W YY=matmul(S,Y)-matmul(zz_1,sigma_T(2:DD+1,t)) XX(:,1:14)=X XX(:,15:28)=matmul(W,X) alpha_T(g,t)=R_g*(A_0**(-1)*alpha_0+(sigma_T(1,t)-sum(sigma_T(2:DD+1,t)*sigma_T(2:DD+1,t)))**(-1.0)*sum((YY-matmul(XX,beta_T(:,t)))))+random_normal()*sqrt(R_g) DEALLOCATE(W); DEALLOCATE(Y) DEALLOCATE(X); DEALLOCATE(XX) DEALLOCATE(S); DEALLOCATE(YY) DEALLOCATE(zz_1) END DO !$omp end parallel !******** Calculate Log likelihood function value *********! Call OMP_SET_NUM_THREADS(Nthreads) !$omp parallel default(shared) private(g,NN,V,C1,C2,C3,C4,W,Y,X,zz_1,zz_2,zero,acc_3v,k,d,i,j, & !$omp & rand,psi_1,psi_2,V1,XX,S1,IPIV,INFO,q_1,q_2,S2,ep_1,ep_2, overage, & !$omp & ep_1v,ep_2v,V2,like_Y1,like_Y2,pp, & !$omp & ithread, iblock, ibegin, iend) ithread = OMP_GET_THREAD_NUM() iblock = (GG+Nthreads-1)/Nthreads ibegin = ithread*iblock+1 iend = min((ithread+1)*iblock,GG) DO g = ibegin,iend NN=N(g) ALLOCATE(V(NN,NN)); ALLOCATE(C1(NN,NN)); ALLOCATE(C2(NN,NN)) ALLOCATE(C3(NN,NN)); ALLOCATE(C4(NN,NN)) ALLOCATE(W(NN,NN)); ALLOCATE(Y(NN)) ALLOCATE(X(NN,14)); ALLOCATE(XX(NN,28)) ALLOCATE(S1(NN,NN)); ALLOCATE(S2(NN,NN)) ALLOCATE(zz_1(NN,DD)); ALLOCATE(zz_2(NN,DD)) ALLOCATE(V1(NN,NN)); ALLOCATE(V2(NN,NN)) ALLOCATE(ep_1(NN)); ALLOCATE(ep_2(NN)) ALLOCATE(ep_1v(NN)); ALLOCATE(ep_2v(NN)) ALLOCATE(IPIV(NN)); ALLOCATE(overage(NN)) ALLOCATE(psi_1(NN,NN)) V=(sigma_T(1,t)-sum(sigma_T(2:DD+1,t)*sigma_T(2:DD+1,t)))*EYE(NN) C1=age(1:NN,1:NN,g) C2=sex(1:NN,1:NN,g) C3=race(1:NN,1:NN,g) C4=sport(1:NN,1:NN,g) W=network(1:NN,1:NN,g) Y=group(1:NN,18,g) X=group(1:NN,1:14,g) overage=group(1:NN,19,g) zz_1=zz(1:NN,1:DD,g) XX(:,1:14)=X XX(:,15:28)=matmul(W,X) IF (missing_num(g)>=1) Y(1:missing_num(g))=Y_T(1:missing_num(g),g) pp=0.0 DO i=1,NN DO j=1,NN IF (j/=i) THEN IF (DD==1) THEN psi_1(i,j)=gamma_T(1,t)+gamma_T(2,t)*c1(i,j)+gamma_T(3,t)*c2(i,j)+gamma_T(4,t)*c3(i,j)+gamma_T(5,t)*overage(i)+gamma_T(6,t)*overage(j) & & -gamma_T(7,t)*abs(zz_1(i,1)-zz_1(j,1)) ELSEIF (DD==2) THEN psi_1(i,j)=gamma_T(1,t)+gamma_T(2,t)*c1(i,j)+gamma_T(3,t)*c2(i,j)+gamma_T(4,t)*c3(i,j)+gamma_T(5,t)*overage(i)+gamma_T(6,t)*overage(j) & & -gamma_T(7,t)*abs(zz_1(i,1)-zz_1(j,1))-gamma_T(8,t)*abs(zz_1(i,2)-zz_1(j,2)) ELSEIF (DD==3) THEN psi_1(i,j)=gamma_T(1,t)+gamma_T(2,t)*c1(i,j)+gamma_T(3,t)*c2(i,j)+gamma_T(4,t)*c3(i,j)+gamma_T(5,t)*overage(i)+gamma_T(6,t)*overage(j) & & -gamma_T(7,t)*abs(zz_1(i,1)-zz_1(j,1))-gamma_T(8,t)*abs(zz_1(i,2)-zz_1(j,2)) -gamma_T(9,t)*abs(zz_1(i,3)-zz_1(j,3)) ELSEIF (DD==4) THEN psi_1(i,j)=gamma_T(1,t)+gamma_T(2,t)*c1(i,j)+gamma_T(3,t)*c2(i,j)+gamma_T(4,t)*c3(i,j)+gamma_T(5,t)*overage(i)+gamma_T(6,t)*overage(j) & & -gamma_T(7,t)*abs(zz_1(i,1)-zz_1(j,1))-gamma_T(8,t)*abs(zz_1(i,2)-zz_1(j,2)) -gamma_T(9,t)*abs(zz_1(i,3)-zz_1(j,3)) & & -gamma_T(10,t)*abs(zz_1(i,4)-zz_1(j,4)) ELSEIF (DD==5) THEN psi_1(i,j)=gamma_T(1,t)+gamma_T(2,t)*c1(i,j)+gamma_T(3,t)*c2(i,j)+gamma_T(4,t)*c3(i,j)+gamma_T(5,t)*overage(i)+gamma_T(6,t)*overage(j) & & -gamma_T(7,t)*abs(zz_1(i,1)-zz_1(j,1))-gamma_T(8,t)*abs(zz_1(i,2)-zz_1(j,2)) -gamma_T(9,t)*abs(zz_1(i,3)-zz_1(j,3)) & & -gamma_T(10,t)*abs(zz_1(i,4)-zz_1(j,4)) -gamma_T(11,t)*abs(zz_1(i,5)-zz_1(j,5)) ENDIF q_1=psi_1(i,j)*W(i,j)-log(1+exp(psi_1(i,j))) pp=pp+q_1 ENDIF END DO END DO S1=eye(NN)-lambda_T(t)*W ep_1=matmul(S1,Y)-matmul(XX,beta_T(:,t))-matmul(zz_1,sigma_T(2:DD+1,t))-alpha_T(g,t) ep_1v=ep_1 V1=V CALL SGESV(NN,1,V1,NN,IPIV,ep_1v,NN,INFO) det1=FindDet(S1, NN) like_Y1=-0.5*NN*log(2*pi)-0.5*NN*log(sigma_T(1,t)-sum(sigma_T(2:DD+1,t)*sigma_T(2:DD+1,t)))+log(det1)-dot_product(ep_1, ep_1v)/2 pp_G(g)=pp+like_Y1 DEALLOCATE(V); DEALLOCATE(C1); DEALLOCATE(C2) DEALLOCATE(C3); DEALLOCATE(C4) DEALLOCATE(W); DEALLOCATE(Y) DEALLOCATE(S1); DEALLOCATE(S2) DEALLOCATE(X); DEALLOCATE(XX) DEALLOCATE(zz_1); DEALLOCATE(zz_2) DEALLOCATE(V1); DEALLOCATE(V2) DEALLOCATE(ep_1); DEALLOCATE(ep_2) DEALLOCATE(ep_1v); DEALLOCATE(ep_2v) DEALLOCATE(IPIV); DEALLOCATE(overage) DEALLOCATE(psi_1) END DO !$omp end parallel logp_T(t)=sum(pp_G) call system_clock ( clock_stop, clock_rate, clock_max ) t_round=t/20 IF (t-t_round*20==0) THEN WRITE(*,*) 'T=', t WRITE(*,'(a6, f12.4)') 'Time:', real ( clock_stop - clock_start ) / real ( clock_rate ) WRITE(*,'(a6,5f12.4)') 'GAMMA:', gamma_T(1:5,t) WRITE(*,'(a6,5f12.4)') 'GAMMA:', gamma_T(6:10,t) WRITE(*,'(a6, f12.4)') 'LAMBDA:', lambda_T(t) WRITE(*,'(a6,6f12.4)') 'SIGMA:', sigma_T(:,t) WRITE(*,'(a6, f12.4)') 'loglike:', logp_T(t) END IF IF (DD==1) THEN WRITE(1000,'(7f12.4)') gamma_T(:,t) WRITE(4000,'(2f12.4)') sigma_T(:,t) ELSEIF (DD==2) THEN WRITE(1000,'(8f12.4)') gamma_T(:,t) WRITE(4000,'(3f12.4)') sigma_T(:,t) ELSEIF (DD==3) THEN WRITE(1000,'(9f12.4)') gamma_T(:,t) WRITE(4000,'(4f12.4)') sigma_T(:,t) ELSEIF (DD==4) THEN WRITE(1000,'(10f12.4)') gamma_T(:,t) WRITE(4000,'(5f12.4)') sigma_T(:,t) ELSEIF (DD==5) THEN WRITE(1000,'(11f12.4)') gamma_T(:,t) WRITE(4000,'(6f12.4)') sigma_T(:,t) ENDIF WRITE(2000,'(f12.4)') lambda_T(t) WRITE(3000,'(28f12.4)') beta_T(:,t) WRITE(6000,'(73f12.4)') alpha_T(:,t) WRITE(7000,'(1f12.4)') logp_T(t) t_round=t/THIN IF (t-t_round*THIN==0) THEN WRITE(5000,'(50f12.4)') zz_T(:,:,:,t/THIN) END IF END DO END PROGRAM SCSAR