Simulations results: good, bad, ugly

[jeremylhour]

Configuration: nb_simu: 1000 # Nb. of simulations sample_size: 10000 # Sample size lambda_x: .9 # Parameter of exponential distribution for X lambda_z: 1 # Parameter of exponential distribution for Z alpha_y: 5 # Parameter of Pareto distribution for Y

b_2+d_2 = .3

Résultats: Theta_0: 2.22

bias:

• smoothed: 0.0001
• standard: 0.0002

MAE:

• smoothed: 0.0037
• standard: 0.0038

RMSE:

• smoothed: 0.0047
• standard: 0.0047

Coverage rate:

• smoothed: 0.8960
• standard: 0.9490

[jeremylhour]

Par contre avec une taille d'échantillon de 1000 aucun problème de taux de couverture...

Configuration: nb_simu: 1000 # Nb. of simulations sample_size: 1000 # Sample size lambda_x: .9 # Parameter of exponential distribution for X lambda_z: 1 # Parameter of exponential distribution for Z alpha_y: 5 # Parameter of Pareto distribution for Y

b_2+d_2 = .3

Résultats: Theta_0: 2.22

bias:

• smoothed: -0.0009
• standard: -0.0001

MAE:

• smoothed: 0.0120
• standard: 0.0121

RMSE:

• smoothed: 0.0149
• standard: 0.0152

Coverage rate:

• smoothed: 0.9500
• standard: 0.9480

[jeremylhour]

Configuration: nb_simu: 1000 # Nb. of simulations sample_size: 100 # Sample size lambda_x: .9 # Parameter of exponential distribution for X lambda_z: 1 # Parameter of exponential distribution for Z alpha_y: 5 # Parameter of Pareto distribution for Y

b_2+d_2 = .3

Résultats: Theta_0: 2.22

bias:

• smoothed: -0.0018
• standard: 0.0027

MAE:

• smoothed: 0.0365
• standard: 0.0379

RMSE:

• smoothed: 0.0462
• standard: 0.0482

Coverage rate:

• smoothed: 0.9270
• standard: 0.9200

[jeremylhour]

Impression

• si b_2+d_2 < 0.5: aucun problème
• si b_2+d_2 > 0.5:
  • si c'est b_2 qui est "grand" alors aucun problème, le taux de couverture > 0.95
  • si c'est d_2 qui est "grand" alors de sérieux problèmes, avec un taux de couverture < 0.95

[jeremylhour]

EDIT: LES RESULTATS NE SONT PAS CORRECTS (PROBLEME DE PARAMETRISATION DE LA LOI EXPONENTIELLE).

Config used:

nb_simu: 10000 # Nb. of simulations sample_size: [100, 200, 500, 1000] # Sample size, can be an array of multiple values lambda_x: [.4, .5, .6, .7, .8] # Parameter of exponential distribution for X lambda_z: [1] # Parameter of exponential distribution for Z alpha_y: [2, 3, 4, 5, 8] # Parameter of Pareto distribution for Y

Runtime: about 2h on SSPdatacould

CICAsymptotics_HugeMCTable.pdf

Commentaires de Xavier:

• comme on l'avait déjà remarqué, il ne semble pas y avoir de grosses différences entre la version lissée et la version non lissée. Ca inciterait plutôt à avoir des résultats sur la version non lissée, si c'est réalisable bien sûr ;
• la vitesse de convergence semble toujours être racine(n). Donc on pourrait conjecturer que quelles que soient les valeurs des paramètres, racinen(n)(theta hat - theta) doit toujours converger vers une certaine distribution. Mais dans certains cas, liés à b_k et d_k, la distribution serait non-normale.
• On a l'impression qu'il n'y a pas beaucoup de cas où on obtient des taux de couverture inférieur à 95%. Du coup on peut se demander si on ne pourrait pas avoir un résultat montrant qu'on pourrait avoir de l'inférence conservatrice sous des hypothèses assez faibles (+ faibles que b_k + d_k < 1/2).

[jeremylhour]

Nouveau résultats:

Config used:

nb_simu: 10000 # Nb. of simulations sample_size: [100, 200, 500, 1000] # Sample size, can be an array of multiple values lambda_x: [.2, .3, .5, .8, .9] # Parameter of exponential distribution for X lambda_z: [1] # Parameter of exponential distribution for Z alpha_y: [1.5, 2, 3, 4, 6, 7, 10] # Parameter of Pareto distribution for Y

CICAsymptotics_HugeMCTable.pdf

Commentaires de Xavier: ce sont des bonnes nouvelles, tout ça. Tes résultats sont peut-être moins "excitants" comme tu dis, mais en effet on les comprend beaucoup mieux !

• Sur ton 2ème point, tu utilises 2 estimateurs de variance différents pour les 2 estimateurs ? Je pose la question car en effet, les 2 estimateurs de theta0 sont assez équivalents en termes de biais et variance donc on ne s'attend pas, si on utilisait le même estimateur de variance, à avoir des taux de couverture très différents.
• Sur ton 3ème point, n=1000 n'est pas si grand, les taux de couverture ne sont quand même pas loin de 95%, j'aurais tendance à attribuer ça en particulier aux diverses approximations asymptotiques, dont l'imperfection de l'estimateur de variance.
• D'ailleurs même si tu ne t'en sors pas dans les intégrales, tu peux estimer précisément la vraie variance asymptotique en tirant un unique mais très très gros échantillon sur lequel tu utilises l'estimateur que tu as programmé.

[jeremylhour]

nb_simu: 10000 # Nb. of simulations sample_size: [100, 500, 1000, 2000] # Sample size, can be an array of multiple values lambda_x: [.2, .3, .5, .8, .9] # Parameter of exponential distribution for X lambda_z: [1] # Parameter of exponential distribution for Z alpha_y: [1.5, 2, 3, 4, 6, 7, 10] # Parameter of Pareto distribution for Y

CICAsymptotics_HugeMCTable.pdf

[jeremylhour]

Même chose que précédemment mais avec l'estimateur de la variance calculé à la Lewbel Schennach

CIC_asymptotics.pdf

[jeremylhour]

nb_simu: 10000 # Nb. of simulations sample_size: [100, 500, 1000, 2000] # Sample size, can be an array of multiple values lambda_x: [.2, .3, .5, .8, .9] # Parameter of exponential distribution for X lambda_z: [1] # Parameter of exponential distribution for Z alpha_y: [1.5, 2, 3, 4, 6, 7, 10] # Parameter of Pareto distribution for Y

Estimateurs:

• 'standard_kernel' : estimateur standard de la FDR de Z, utilisation d'un estimateur par noyau pour f_Y(.)
• 'standardxavier': estimateur standard de la FDR de Z, utilisation de la technique proposée par Xavier (utiliser X(i+1) et X(i-1) pour la dérivée en X(i)) pour estimer 1/f_Y(F_Y(.))
• 'smooth_kernel': estimateur smooth de la FDR de Z, utilisation d'un estimateur par noyau pour f_Y(.)
• 'smoothls': estimateur smooth de la FDR de Z, adaptation de la méthode de Lewbel et Schennach (utiliser X(i+1) et X(i) pour la dérivée en (X(i+1)+X_(i))/2 et suppression des dupliqués pour éviter la division par zéro) pour estimer 1/f_Y(F_Y(.))
• 'smoothxavier': estimateur smooth de la FDR de Z, utilisation de la technique proposée par Xavier (utiliser X(i+1) et X(i-1) pour la dérivée en X(i)) pour estimer 1/f_Y(F_Y(.))

Rappel: la FDR de Z est utilisé pour le calcul de \hat U = \hat F_Z(X). L'estimateur de f_Y(.) n'entre en comtpe que dans le calcul de la variance.

Attention: le code a mis 1,5 jours à tourner environ.

CICAsymptotics_HugeMCTable.pdf

[jeremylhour]

Nouveaux résultats avec deux DGP, et en affichant la taille des intervalles de confiance. J'ai réduit le nombre de tirage à 5000 et je n'ai pas fait pour des tailles d'échantillon égales à 2000 pour limiter le temps de calcul.

1) DGP "habituel". Ce sont les mêmes résultats que précédemment, simplement on peut y voir la taille des IC. On voit que "smooth_kernel" a un gros problème. Sur tous les DGP.

2) Le nouveau DGP avec des lois normales. On a:

• Y ~ N(0,1),
• Z ~ N(0,1),
• X ~ N(mu_x, variance_x). Donc à chaque bloc, c'est mu_x ou variance_x qui varie. D'après mes calculs variance_x > variance_z ne pose pas vraiment de problème car la densité f_U tend vers zéro sur les bords -- je n'ai fait aucune simulation pour ce cas pour le moment. En revanche, j'avais estimé que variance_z > variance_x pouvait poser problème dans certains cas selon que variance_x/variance_z soit plus ou moins grand que 1/2 (d'où les bornes pages 6 de notre document). J'ai l'impression que ça ne pose pas plus problème que ça. En revanche, l'estimation est très mauvaise quand mu_x est très grand par rapport à mu_z. Or, d'après mes calculs mu_x-mu_z ne devant plus jouer aux bords (U proche de 0 ou de 1). Est-ce que c'est parce que du coup ça envoit les X dans la queue de distribution des Z, et comme il n'y a pas beaucoup d'observations dans cette zone ça biaise beaucoup les résultats? Ce serait mon intuition "à chaud".

CICAsymptotics_HugeMCTable.pdf

[jeremylhour]

There was an overflow with estimator 'smooth_kernel' especially when considering the Exponential DGP. This is not a bug. The standard error is very large, and it comes from the 'inverse of the density' terms that goes into the standard error. With that DGP, some points are very extreme, so the density estimator becomes very close to zero at the right tail, which becomes a problem when a \hat U is large and falls, e.g., between two extreme values.

Code points for the problem :

• EXPONENTIAL_simulations.py, line 110 :

  theta_smooth, sigma_smooth = estimator_unknown_ranks(y, x, z, method="smoothed", se_method="kernel")

• func_main.py, line 68 :

  if se_method == "kernel":
  inv_density = 1/kernel_density_estimator(x=np.quantile(y, u_hat), data=y)

[jeremylhour]

We tried the following config for the Gaussian DGP :

nb_simu: 100 # Nb. of simulations sample_size: [1000, 10000, 20000] # Sample size, can be an array of multiple values mu_x: [0, 3, 5] # mean of X (Gaussian distribution) variance_x: [.5] # Variance of X (Gaussian distribution)

When mu_X = 5 and sample_size is 20000, the code crashes as there seem to be a lot of overflow. Results are reported in the document for the rest. For mu_X = 3, results on coverage rate seems to go in the right direction as the sample size increases. Although it seems that to get close to .95, a lot of data is required. We also tried with a sample size of 100,000 but even one iteration would not finish in a day.

DGP_Gaussian_Large_Sample.pdf

[jeremylhour]

[jeremylhour]