subroutine statss(ns,k,a,b,c,alp,nwant,dwant,lrand,randstat,nper,nsim,lstat, & stat) ! ! Compute statistics for the estimated model: ! ! S[t+1] = A S[t] + B eps[t+1] , Eeps eps'=I ! Y[t] = C S[t] ! ! Statistics reported: ! ! (1) A0 from a SVAR of Y[t] regressed on infinite lags, ! ! Y[t] = B1 Y[t-1] + M B1 Y[t-2] + M^2 B1 Y[t-3]+ ... + e[t] ! = B1 Y[t-1] + B2 Y[t-2] + B3 Y[t-3]+ ... + e[t] ! ! A0 = C*B ! ! ==> stat(1:k*k) = [A0(1,:),A0(2,:),...,A0(k,:)] ! ! (2) Population variance decomposition for Y ! ! ==> stat(k*k+1:2k*k+k) = [var(Y[t]|all shocks) ! var(Y[t]|shock 1 only)/var(Y[t]|all) ! : ! var(Y[t]|shock k only)/var(Y[t]|all) ! ! (3) Second moments for variables of interest, W[t], where ! ! W[t] = D* X[t], D is n x k ! Y[i,t] = (1-alp(i)L) X[i,t] ! ! and variables in W[t] are HP-filtered. The statistics reported ! are standard deviations, autocorrelations, contemporaneous ! correlations, and the variance decomposition of HP-filtered W[t] ! for nsim simulations of length nper periods for the model. ! ! ==> stat((2k+1)*k+1:lstat) = [std(HP W[t]) ! cor(HP W[t], HP W[t-1]) ! cor(HP W[:,t], HP W[1,t]) ! : ! cor(HP W[:,t], HP W[nwant,t])] ! var(HP W[t]|shock 1 only)/var(.|all) ! : ! var(HP W[t]|shock k only)/var(.|all)] ! ! where lstat = (2k+1)*k+(nwant+k+2)*nwant ! implicit none integer, parameter :: lwork=5000 integer, intent(in) :: ns,k,nwant,lrand,nper,nsim,lstat integer :: i,j,l,info,n,s,n1,n2 integer, dimension(lwork) :: iwork real, dimension(ns,ns), intent(in) :: a real, dimension(ns,k), intent(in) :: b real, dimension(k,ns), intent(in) :: c real, dimension(lrand), intent(in) :: randstat real, dimension(lstat), intent(out) :: stat real, dimension(k,k) :: a0,vy,vj,epsi real, dimension(ns,k) :: tem1 real, dimension(ns,ns) :: at,vs,atem,tem2 real, dimension(ns,1) :: state real, dimension(k,1) :: y,eps,x,alp real, dimension(nwant,1) :: wtem real, dimension(nwant) :: sum1,sum2,sum3,sum4,sum5,sum6, & sum7,mean,meanf,meanl,std,var, & varf,varl,acov,acor real, dimension(nsim,nwant) :: varhld real, dimension(nwant*nwant) :: sum8,cor real, dimension(nwant,nper) :: w real, dimension(nwant,k) :: dwant real, dimension(nper) :: wn,wt,wd real, dimension((nwant+2+k)*nwant) :: kp real, dimension(lwork) :: work stat = 0. if (lstat<(2*k+1)*k+(nwant+2)*nwant+nwant*k) then write(*,*) 'ERROR: In callss.f90, fix nstat in call to stats.f90' endif !======================== ! First set of statistics !======================== ! ! A0 = C*B ! a0 = matmul(c,b) l = 1 do i=1,k do j=1,k stat(l) = a0(i,j) l = l+1 enddo enddo !========================= ! Second set of statistics !========================= ! ! Compute Vs and Vy that satisfy Vs = A*Vs*A' + B*B' ! Vy = C*Vs*C' ! vs = matmul(b,transpose(b)) at = -transpose(a) atem = a call sb04qd(ns,ns,atem,ns,at,ns,vs,ns,tem2,ns,iwork,work,lwork,info) tem1 = matmul(vs,transpose(c)) vy = matmul(c,tem1) l = k*k+1 do i=1,k stat(l) = vy(i,i) l = l+1 enddo ! ! Variance decomposition (V(i) due to each eps(i)) ! V(i) = A*V(i)*A' + B*epsi*B', epsi has 1 in i,i element ! l = k*k+k+1 do j=1,k epsi = 0. epsi(j,j) = 1. tem1 = matmul(b,epsi) vs = matmul(tem1,transpose(b)) atem = a at = -transpose(a) call sb04qd(ns,ns,atem,ns,at,ns,vs,ns,tem2,ns,iwork,work,lwork,info) tem1 = matmul(vs,transpose(c)) vj = matmul(c,tem1) do i=1,k stat(l) = vj(i,i)/vy(i,i) l = l+1 enddo enddo !========================= ! Third set of statistics !========================= ! ! KP Statistics for the variables of interest (W) ! ! Example: Suppse Y and W are given by ! ! Y[t] = [ log (y[t]/l[t]) - log(y[t-1]/l[t-1]) ] ! [ log l[t] - alpha log l[t-1] ] ! ! W[t] = [ log y[t] ] ! [ log l[t] ] ! ! Then: alp = [1,alpha]' ! D = [1 1] ! [0 1] ! nwant = number of rows of W ! = 2 ! state = 0. kp = 0. l = 1 n1 = (nwant+2)*nwant do s=1,nsim x = 0. do j=1,nper do i=1,k eps(i,1) = randstat(l) l = l+1 enddo state = matmul(a,state)+matmul(b,eps) y = matmul(c,state) x = alp*x -y wtem = matmul(dwant,x) w(:,j) = wtem(:,1) enddo do j=1,nwant wn = w(j,:) call hptrend(nper,wn,1600., wt,wd) w(j,:) = wd enddo ! ! With HP-filtered data, compute second moments (as in Kydland-Prescott) ! sum1 = 0. sum2 = 0. sum3 = 0. sum4 = 0. sum5 = 0. sum6 = 0. sum7 = 0. sum8 = 0. do j=1,nper sum1 = sum1 + w(:,j) sum2 = sum2 + w(:,j)*w(:,j) enddo do j=2,nper sum3 = sum3 + w(:,j) sum4 = sum4 + w(:,j-1) sum5 = sum5 + w(:,j)*w(:,j-1) sum6 = sum6 + w(:,j)*w(:,j) sum7 = sum7 + w(:,j-1)*w(:,j-1) enddo mean = sum1/float(nper) meanf = sum3/float(nper-1) meanl = sum4/float(nper-1) var = (sum2-float(nper)*mean*mean)/float(nper-1) std = sqrt(var) acov = (sum5-float(nper-1)*meanf*meanl)/float(nper-2) varf = (sum6-float(nper-1)*meanf*meanf)/float(nper-2) varl = (sum7-float(nper-1)*meanl*meanl)/float(nper-2) acor = acov/(sqrt(varf)*sqrt(varl)) do i=1,nwant do n=1,nwant do j=1,nper sum8((i-1)*nwant+n) = sum8((i-1)*nwant+n)+w(i,j)*w(n,j) enddo cor((i-1)*nwant+n) = (sum8((i-1)*nwant+n)-float(nper)*mean(i)* & mean(n))/(float(nper-1)*std(i)*std(n)) enddo enddo kp(1:nwant) = kp(1:nwant)+std kp(nwant+1:2*nwant) = kp(nwant+1:2*nwant)+acor kp(2*nwant+1:n1) = kp(2*nwant+1:n1)+cor varhld(s,:) = var enddo do n=1,k state = 0. l = 1 do s=1,nsim x = 0. do j=1,nper do i=1,k if (i==n) then eps(i,1) = randstat(l) else eps(i,1) = 0. endif l = l+1 enddo state = matmul(a,state)+matmul(b,eps) y = matmul(c,state) x = alp*x -y wtem = matmul(dwant,x) w(:,j) = wtem(:,1) enddo do j=1,nwant wn = w(j,:) call hptrend(nper,wn,1600., wt,wd) w(j,:) = wd enddo ! ! With HP-filtered data, compute second moments (as in Kydland-Prescott) ! sum1 = 0. sum2 = 0. do j=1,nper sum1 = sum1 + w(:,j) sum2 = sum2 + w(:,j)*w(:,j) enddo mean = sum1/float(nper) var = (sum2-float(nper)*mean*mean)/float(nper-1) n1 = (nwant+2)*nwant+(n-1)*nwant+1 n2 = (nwant+2)*nwant+n*nwant kp(n1:n2) = kp(n1:n2)+var/varhld(s,:) enddo enddo stat((2*k+1)*k+1:lstat) = kp/float(nsim) end subroutine statss