📄 gmm.f90
字号:
write(10,*)W(i,j)*real(nobs,kind(0.0d0))
end do
end do
close(10)
as(1,:)=a
call residual(as,resid,nerr)
open(10,file='series.out',status='old')
rewind(10)
do i=1,no_of_equations
write(10,*)'equation',i
do j=1,nobs
write(10,*)indep_var(i,j),resid(1,i,j)
end do
end do
close(10)
RETURN
CONTAINS
SUBROUTINE WRITE_STATS(a,gmm_stage)
! WRITE ASYMPTOTIC COVARIANCE MATRIX (transpose(D)WD)^-1
logical gmm_stage
real(8), dimension(:) ::a
real(8), dimension(1:1,size(a)) ::as
real(8), dimension(max_num_moments) ::Mean,second_M, &
third_M,fourth_M
real(8) tmp0,tmp1,tmp2,tmp3,tmp4
CALL calccovmtx(a,deriv_type,eps_vector,covmtx,num_moments, &
matrix_ill)
if(matrix_ill) then
write(*,*)'inverse of covariance matrix is ill conditioned'
write(*,*)'statistics are skipped'
matrix_ill=.false.
goto 200
end if
write(*,*)'covariance matrix for parameters:'
write(*,*) covmtx
write(*,*)
write(*,*)'parameters and their t-statistics'
do i=1,size(a)
write(*,*)a(i),a(i)/sqrt(covmtx(i,i))
end do
200 as(1,:)=a
call residual(as,resid,nerr)
!! Calculate some basic stats
write(*,*)
do i=1,no_of_equations
write(*,*)'Basic statistics for equation................',i
tmp0=sum(indep_var(i,:))/real(nobs,kind(0.0d0))
write(*,*)'Mean of independent variable ',tmp0
tmp1=(dot_product(indep_var(i,:),indep_var(i,:))/ &
real(nobs,kind(0.0d0)))-tmp0*tmp0
write(*,*)'Var(indep. variable) around mean:divisor N',tmp1
tmp2=sum(resid(1,i,:))/real(nobs,kind(0.0d0))
write(*,*)'Mean of residual ',tmp2
tmp3=(dot_product(resid(1,i,:),resid(1,i,:))/ &
real(nobs,kind(0.0d0)))
write(*,*)'Var(residual) around zero :divisor N ',tmp3
tmp4=tmp3-tmp2*tmp2
write(*,*)'Var(residual) around mean :divisor N ',tmp4
tmp4=1.0d0 - tmp3/tmp1
write(*,*)'R-squared statistic ',tmp4
tmp0=1.0d0-(1.0d0-tmp4)*real(nobs-1,kind(0.0d0))/ &
real(nobs-np,kind(0.0d0))
write(*,*)'adjusted R-squared statistic ',tmp0
end do
write(*,*)
if(gmm_stage) then
write(*,*)
write(*,*)'Test of overidentifying restrictions:'
write(*,*)'J-statistic = ',fret
write(*,*)'distributed as Chi squared with'
write(*,*) num_moments-size(a),' D.F.'
write(*,*)'Weighting matrix type = ',weighting_type
write(*,*)'with nma =',nma
as(1,:)=a
! the sample moment vector is the average of the row
! vectors in FMTX
call CONSTFMTX(as)
mean(1:num_moments)=sum(Fmtx(1,1:num_moments, &
:),dim=2)/real(nobs,kind(0.0d0))
Fmtx(1,1:num_moments,:)=Fmtx(1,1:num_moments,:)-spread( &
Mean(1:num_moments),dim=2,ncopies=nobs)
second_M(1:num_moments) = sum(Fmtx(1,1:num_moments, &
:)**2,dim=2)/real(nobs-1,kind(0.0d0))
third_M(1:num_moments) = sum(Fmtx(1,1:num_moments, &
:)**3,dim=2)/real(nobs-1,kind(0.0d0))
fourth_M(1:num_moments) = sum(Fmtx(1,1:num_moments, &
:)**4,dim=2)/real(nobs-1,kind(0.0d0))
write(*,*)
write(*,*) 'Central moments of orthogonality conditions:'
write(*,*)' MEAN 2nd ', &
' 3rd 4th'
do i=1,num_moments
write(*,*) mean(i),second_M(i),third_M(i),fourth_M(i)
end do
end if
write(*,*)
write(*,*)
write(*,*)
write(*,*)
return
END SUBROUTINE WRITE_STATS
SUBROUTINE cross_equ_cov(a,crosscovmtx)
! subroutine to calculate cross equation covariance after 2SLS
real(8), dimension(:) ::a
real(8), dimension(:,:) ::crosscovmtx
real(8), dimension(1:1,size(a)) ::as
as(1,:)=a
call residual(as,resid,nerr)
if(nerr.lt.0) write(*,*) 'error in residual calculation'
do j=1,no_of_equations
Fmtx(1,j,:)=resid(1,j,:) ! use Fmtx as work space
end do
crosscovmtx=matmul(Fmtx(1,1:no_of_equations,:),transpose( &
Fmtx(1,1:no_of_equations,:)))/real(nobs- &
size(a),kind(0.0d0))
RETURN
END SUBROUTINE cross_equ_cov
END SUBROUTINE GMM
! *********************************************************
SUBROUTINE setinst
! This subroutine sets up the instruments in FMTX
! A row of Fmtx will contain, after this, the observed instruments
! for the corresponding errors; therefore nobs columns in Fmtx.
! There are no_of_equations equations and num_inst(i) instruments
! for equation i, consequently there are 'SUM' rows in Fmtx where
! SUM is:
! do j=1,nobs
! do i=1,no_of_equations
! SUM=SUM+num_inst(i)
! end do
! end do
!
! eg. fmtx= z11i
! z12i
! z13i
! z21i
! z22i
implicit none
include 'dim_dat.inc'
integer, parameter::max_num_moments=no_of_equations* &
no_of_instruments
integer, parameter ::dim1_of_Fmtx=max(max_num_as,np*5)
real(8), dimension(dim1_of_Fmtx,max_num_moments,nobs) ::Fmtx
real(8), dimension(max_num_moments,max_num_moments) ::W
real(8), dimension(no_of_instruments,nobs) ::Z
integer, dimension(no_of_equations) ::num_inst
integer i,m,j
character(len=no_of_instruments),dimension(no_of_equations) &
::instruments
common/instruments/ Z
common/flag_instruments/ instruments
common/gmm_specific/ Fmtx,W,num_inst
m=0
do j=1,no_of_equations
do i=1,no_of_instruments
if(instruments(j)(i:i).eq.'1') then
m=m+1
Fmtx(1,m,:)=Z(i,:)
end if
end do
end do
return
end subroutine setinst
! *********************************************************
SUBROUTINE CONSTFMTX(as)
! This subroutine constructs FMTX.
! A row of Fmtx will contain, after this, the observed instruments
! multiplied by the corresponding error; therefore there are nobs
! columns in Fmtx. There are no_of_equations equations and
! num_inst(i) instruments
! for equation i, consequently there are 'SUM' rows in Fmtx where
! SUM is:
! do j=1,nobs
! do i=1,no_of_equations
! SUM=SUM+num_inst(i)
! end do
! end do
!
! eg. fmtx= e1i*z11i
! e1i*z12i
! e1i*z13i
! e2i*z21i
! e2i*z22i
implicit none
include 'dim_dat.inc'
integer, parameter ::max_num_moments=no_of_equations* &
no_of_instruments
integer, parameter ::dim1_of_Fmtx=max(max_num_as,np*5)
real(8), dimension(dim1_of_Fmtx,max_num_moments,nobs) ::Fmtx
real(8), dimension(max_num_moments,max_num_moments) ::W
integer, dimension(no_of_equations) ::num_inst
real(8), dimension(:,:) ::as
real(8), dimension(size(as,dim=1),no_of_equations,nobs)::resid
integer m,j,nerr,no_as,n
common/gmm_specific/ Fmtx,W,num_inst
interface
SUBROUTINE residual(as,resid,nerr)
real(8), dimension(:,:) ::as
real(8), dimension(:,:,:) ::resid
integer nerr
end SUBROUTINE residual
end interface
! put the instruments into Fmtx
no_as=size(as,dim=1)
CALL setinst
! if(no_as.gt.1) Fmtx(2:no_as,:,:)=spread(Fmtx(1,:,:),dim=1, &
! ncopies=no_as-1)
! replace this command with the following
! due to machine constraints
do m=2,no_as
Fmtx(m,:,:)=Fmtx(1,:,:)
end do
! create the Fmtx by getting a vector of residuals and multiplying
! by the appropriate instruments
m = 1
call residual(as,resid,nerr)
do j=1,no_of_equations
! Fmtx(1:no_as,m:m+num_inst(j)-1,:)=spread(resid(1: &
! no_as,j,:),dim=2,ncopies=num_inst(j))* &
! Fmtx(1:no_as,m:m+num_inst(j)-1,:)
! the above line is equivalent to the following 3 lines
! (use if further machine constraints):
do n = m , (m+num_inst(j)-1)
Fmtx(1:no_as,n,:)=resid(1:no_as,j,:)*Fmtx(1:no_as,n,:)
end do
m=m+num_inst(j)
end do
RETURN
END SUBROUTINE CONSTFMTX
! *********************************************************
SUBROUTINE calccovmtx(a,deriv_type,eps_vector,covmtx, &
num_moments,matrix_ill)
! This subroutine calculates the covariance matrix for the
! parameters (a vector)
!
! Note: at this point in the GMM procedure, the Fmtx is no longer
! required. Consequently the space will be used for some of the
! calculations.
!
implicit none
include 'dim_dat.inc'
integer, parameter ::max_num_moments=no_of_equations* &
no_of_instruments
integer, parameter ::dim1_of_Fmtx=max(max_num_as,np*5)
real(8), dimension(dim1_of_Fmtx,max_num_moments,nobs) ::Fmtx
real(8), dimension(max_num_moments,max_num_moments) ::W
integer, dimension(no_of_equations) ::num_inst
integer num_moments
integer m,j,i,n,deriv_type
real(8), dimension(:) ::a,eps_vector
real(8), dimension(:,:) ::covmtx
real(8), dimension(size(a),no_of_equations,nobs) ::grad
real(8), dimension(max_num_moments,size(a)) ::Dmtx
logical matrix_ill
common/gmm_specific/ Fmtx,W,num_inst
include 'blas.inc'
interface
subroutine GETGRAD(a,grad,deriv_type,eps_vector)
real(8), dimension(:) ::a,eps_vector
real(8), dimension(:,:,:) ::grad
integer deriv_type
end subroutine getgrad
end interface
! put the instruments into Fmtx
CALL setinst
m = 1
call getgrad(a,grad,deriv_type,eps_vector)
do j=1,no_of_equations
do i=1,size(a)
do n=m,(m+num_inst(j)-1)
Dmtx(n,i)=dot_product(grad(i,j,:),Fmtx(1,n,:))
end do
end do
m=m+num_inst(j)
end do
covmtx=matmul(matmul(transpose(Dmtx(1:num_moments,:)), &
W(1:num_moments,1:num_moments)), &
Dmtx(1:num_moments,:))
call INVERSE(covmtx,matrix_ill)
RETURN
END SUBROUTINE calccovmtx
! *********************************************************
SUBROUTINE funcgmm0(spare_params,as,F,nerr)
! Defines the GMM objective function
! It is of the form required by MINIMIZE.
! spare_params = integer parameters which need to be passed.
! (1) = weighting_type
! (2) = nma
! (3) = num_moments
! as = set of different parameter vectors; as(1,:),
! as(2,:) , ..., as(i,:) represents the set of
! parameter vectors which the objective function is
! calculated for.
! F = Objective function value
! nerr = returns a -ve value if an error occurs.
implicit none
include 'dim_dat.inc'
integer, parameter ::max_num_moments=no_of_equations* &
no_of_instruments
integer, parameter ::dim1_of_Fmtx=max(max_num_as,np*5)
real(8), dimension(dim1_of_Fmtx,max_num_moments,nobs) ::Fmtx
real(8), dimension(max_num_moments,max_num_moments) ::W
integer, dimension(no_of_equations) ::num_inst
integer, dimension(:) ::spare_params
integer nerr,no_as,i
real(8), dimension(:,:) ::as
real(8), dimension(size(as,dim=1),spare_params(3)) ::M,MTW
real(8), dimension(:) ::F
common/gmm_specific/ Fmtx,W,num_inst
interface
SUBROUTINE constfmtx(as)
real(8), dimension(:,:) ::as
END SUBROUTINE constfmtx
END INTERFACE
! num_moments=spare_params(3)
! the sample moment vector is the average of the row vectors in FMTX
no_as=size(as,dim=1)
call CONSTFMTX(as)
M(1:no_as,:)=sum(Fmtx(1:no_as,1:spare_params(3),:),dim=3)
! the GMM objective is to minimize: transpose(M).W.M where W
! is a weighting matrix
⌨️ 快捷键说明
复制代码
Ctrl + C
搜索代码
Ctrl + F
全屏模式
F11
切换主题
Ctrl + Shift + D
显示快捷键
?
增大字号
Ctrl + =
减小字号
Ctrl + -