adding tests for the inference part

src/econometrics.jl

@@ -44,7 +44,9 @@ function calculate_Avar!(sMMProblem::SMMProblem, theta0::Array{Float64,1}, tData
 
   # sandwich formula, where tau capture the extra noise introduced by simulation,
   # compared to the usual GMM estimator
-  sMMProblem.Avar = (1+tau)*inv_Sigma1*Sigma2*inv_Sigma1
+  # I use the function Symmetric, because otherwise a test issymmetric()
+  # on the asymptotic variance will return false because of numerical approximations
+  sMMProblem.Avar = Symmetric((1+tau)*inv_Sigma1*Sigma2*inv_Sigma1)
 
 end
 

test/runtests.jl

@@ -837,3 +837,205 @@ end
 
 
 end
+
+
+@testset "Testing Inference" begin
+
+  tolLinear = 0.05
+
+  # Inference in the linear model
+  #------------------------------
+  srand(1234)         #for replicability reasons
+  T = 100000          #number of periods
+  P = 2               #number of dependent variables
+  beta0 = rand(P)     #choose true coefficients by drawing from a uniform distribution on [0,1]
+  alpha0 = rand(1)[]  #intercept
+  theta0 = 0.0        #coefficient to create serial correlation in the error terms
+  println("True intercept = $(alpha0)")
+  println("True coefficient beta0 = $(beta0)")
+  println("Serial correlation coefficient theta0 = $(theta0)")
+
+  # Simulation of a sample:
+  # Generation of error terms
+  #--------------------------
+  # row = individual dimension
+  # column = time dimension
+  U = zeros(T)
+  d = Normal()
+  U[1] = rand(d, 1)[] #first error term
+  # loop over time periods
+  for t = 2:T
+      U[t] = rand(d, 1)[] + theta0*U[t-1]
+  end
+  # Let's simulate x_t
+  #-------------------
+  x = zeros(T, P)
+
+  d = Uniform(0, 5)
+  for p = 1:P
+          x[:,p] = rand(d, T)
+  end
+
+  # Let's calculate the resulting y_t
+  y = zeros(T)
+
+  for t=1:T
+      y[t] = alpha0 + x[t,1]*beta0[1] + x[t,2]*beta0[2] + U[t]
+  end
+
+  myProblem = SMMProblem(options = SMMOptions(maxFuncEvals=500, saveSteps = 500, globalOptimizer = :dxnes, localOptimizer = :LBFGS, minBox = false, showDistance = false));
+
+  # Empirical moments
+  #------------------
+  dictEmpiricalMoments = OrderedDict{String,Array{Float64,1}}()
+  dictEmpiricalMoments["mean"] = [mean(y); mean(y)] #informative on the intercept
+  dictEmpiricalMoments["mean_x1y"] = [mean(x[:,1] .* y); mean(x[:,1] .* y)] #informative on betas
+  dictEmpiricalMoments["mean_x2y"] = [mean(x[:,2] .* y); mean(x[:,2] .* y)] #informative on betas
+  dictEmpiricalMoments["mean_x1y^2"] = [mean((x[:,1] .* y).^2); mean((x[:,1] .* y).^2)] #informative on betas
+  dictEmpiricalMoments["mean_x2y^2"] = [mean((x[:,2] .* y).^2); mean((x[:,2] .* y).^2)] #informative on betas
+
+  set_empirical_moments!(myProblem, dictEmpiricalMoments)
+
+  dictPriors = OrderedDict{String,Array{Float64,1}}()
+  dictPriors["alpha"] = [0.5, 0.001, 1.0]
+  dictPriors["beta1"] = [0.5, 0.001, 1.0]
+  dictPriors["beta2"] = [0.5, 0.001, 1.0]
+
+  set_priors!(myProblem, dictPriors)
+
+  # x[1] corresponds to the intercept
+  # x[1] corresponds to beta1
+  # x[3] corresponds to beta2
+  @everywhere function functionLinearModel(x; nbDraws::Int64 = 1000000, burnInPerc::Int64 = 10)
+
+      # Structural Model
+      #-----------------
+      srand(1234) #for replicability reasons
+      T = nbDraws
+      P = 2       #number of dependent variables
+
+      alpha = x[1]
+      beta = x[2:end]
+      theta = 0.0     #coefficient to create serial correlation in the error terms
+
+
+      # Creation of error terms
+      # row = individual dimension
+      # column = time dimension
+      U = zeros(T)
+      d = Normal()
+      U[1] = rand(d, 1)[] #first error term
+      # loop over time periods
+      for t = 2:T
+          U[t] = rand(d, 1)[] + theta*U[t-1]
+      end
+
+      simX = zeros(T, P)
+      d = Uniform(0, 5)
+      for p = 1:P
+              simX[:,p] = rand(d, T)
+      end
+
+      # Let's calculate the resulting y_t
+      y = zeros(T)
+
+      for t=1:T
+          y[t] = alpha + simX[t,1]*beta[1] + simX[t,2]*beta[2] + U[t]
+      end
+
+      # Get rid of the burn-in phase:
+      #------------------------------
+      startT = div(nbDraws, burnInPerc)
+
+      # Moments:
+      #---------
+      output = OrderedDict{String,Float64}()
+      output["mean"] = mean(y[startT:nbDraws])
+      output["mean_x1y"] = mean(simX[startT:nbDraws,1] .* y[startT:nbDraws])
+      output["mean_x2y"] = mean(simX[startT:nbDraws,2] .* y[startT:nbDraws])
+      output["mean_x1y^2"] = mean((simX[startT:nbDraws,1] .* y[startT:nbDraws]).^2)
+      output["mean_x2y^2"] = mean((simX[startT:nbDraws,2] .* y[startT:nbDraws]).^2)
+
+      return output
+  end
+
+  set_simulate_empirical_moments!(myProblem, functionLinearModel)
+
+  # Construct the objective function using:
+  #* the function: parameter -> simulated moments
+  #* emprical moments values
+  #* emprical moments weights
+  construct_objective_function!(myProblem)
+
+  # Run the optimization in parallel using n different starting values
+  # where n is equal to the number of available workers
+  #--------------------------------------------------------------------
+  @time listOptimResults = local_to_global!(myProblem, verbose = true)
+
+  minimizer = smm_local_minimizer(myProblem)
+
+  # The minimizer should not be too far from the true values
+  #---------------------------------------------------------
+  @test minimizer[1] ≈  alpha0[1] atol = tolLinear
+
+  @test minimizer[2] ≈  beta0[1] atol = tolLinear
+
+  @test minimizer[3] ≈  beta0[2] atol = tolLinear
+
+  # Empirical Distance matrix
+  #--------------------------
+  X = zeros(T, 5)
+
+  X[:,1] = y
+  X[:,2] = (x[:,1] .* y)
+  X[:,3] = (x[:,2] .* y)
+  X[:,4] = (x[:,1] .* y).^2
+  X[:,5] = (x[:,2] .* y).^2
+
+  Sigma0 = cov(X)
+
+  set_Sigma0!(myProblem, Sigma0)
+
+  @test myProblem.Sigma0 == Sigma0
+
+  nbDraws = 1000000 #number of draws in the simulated data
+  calculate_Avar!(myProblem, minimizer, T, nbDraws)
+
+  # The asymptotic variance should be
+  # * symmetric
+  # * positive semi-definite
+  #-----------------------------------
+  @test issymmetric(myProblem.Avar) == true
+
+  # test for positive semi-definiteness
+  # all eigenvalues should be non-negative
+  eigv = eigvals(myProblem.Avar)
+
+  counterNegativeEigVals = 0
+  for i=1:length(eigv)
+    if eigv[i] < 0
+      counterNegativeEigVals += 1
+    end
+  end
+
+  @test counterNegativeEigVals == 0
+
+  # summary table:
+  #---------------
+  df = summary_table(myProblem, minimizer, T, 0.05)
+
+  # first column : point estimates
+  @test df[:Estimate] == minimizer
+
+  # 2nd column : std error
+  for i =1:size(df,1)
+    @test df[:StdError][i] > 0.
+  end
+
+  # confidence interval
+  for i =1:size(df,1)
+    @test df[:ConfIntervalLower][i] <= df[:ConfIntervalUpper][i]
+  end
+
+
+end