Implementação em C# da Distribuição Normal

A distribuição normal (ou Gaussiana) tende a entrar em tudo, mesmo, conforme Nassin Taleb, onde não deveria. É uma distribuição simples, facilmente integrável e derivável, e aparece naturalmente do teorema do limite central.

Os praticantes de matemática financeira frequentemente precisam computar a Função Cumulativa ou a Inversa da Cumulativa da Distribuição normal. Por exemplo no cálculo de fatores de confiança paramétrico para o cálculo de alguns VaR Matriciais quanto na estimativa do prêmio de opções, que, via de regra, assumem distribuição normal (ou Log-Normal) de suas variáveis. Não existem fórmulas fechadas para o cálculo destas funções, então aproximações numéricas são necessárias.

No Excel usa-se NormSDist e NormSInv, (Dist.NormP e Inv.NormP no Excel em Português). NormSInv, em especial tem má fama devido a implementação que existia até o Excel 2000, que era de precisão simples (“8 digitos”, single) apesar do Excel tratar todo número como duplo (“16 digitos”, double).

Isso mudou com o Excel 2003, quando a Microsof, sob crítica, finalmente trocou o algoritmo para um de dupla precisão.

Artigo muito citado de  Graeme West ilustra esses problemas a partir de algumas criticas que Espen Haug recebeu por distribuir seu famoso livro,  The Complete Guide to Option Pricing Formulas, com uma implementação simples da Cumulativa chamada Abramowitz & Stegun 1974, que funciona bem, mas tem precisão simples, apenas 8 digitos. Essa implementação de pouca precisão, quando aplicada em algumas das formulações de opções, levava a premios negativos, um óbvio absurdo.

No Brasil, em opções sobre IDI, que são Vanilla, a aplicação de uma formulação de pouca precisão leva a distorções perceptíveis no preço, dado o montante do Spot e do Strike, hoje da ordem de 300.000.

Graeme sugere usar um algoritmo de precisão double, Hart 1968, para sanar os problemas que uma má implementação de NormSDist pode causar.

E para calcular a função inversa, NormSInv? Pode-usar um algoritmo de localização (como Newton, Bissecção, etc) para achar numericamente o valor, e é assim que o Excel faz. No entanto isso vai requerer mais poder computacional, o que pode tornar-se critico ao usar-se NormSInv, por exemplo, numa simulação Monte-Carlo. Uma aproximação numérica, de dupla precisão, pode ser encontrada em http://lib.stat.cmu.edu/apstat/241, mas em Fortran.

Movido por estas questões contruímos uma pequena classe utilitária para o cálculo destas funções mais comuns, traduzindo os algoritmo de Hart 1968 e do AS241 para C#. Nos foi útil, pode ser para outros.

Distribuição Normal

1. /// <summary>
2.     /// Funções da distribuição normal standard
3.     /// </summary>
4.     public static class StandardNormalDistribution
5.     {
6.         /// <summary>
7.         /// Alias para o método StandardCumulative desta classe, com o mesmo nome que a função equivalente do Excel em inglês
8.         /// </summary>
9.         /// <param name="z">Is the value for which you want the distribution</param>
10.         /// <returns>
11.         /// Returns the standard normal cumulative distribution function.
12.         /// The distribution has a mean of 0 (zero) and a standard deviation of one.
13.         /// Use this function in place of a table of standard normal curve areas.
14.         /// </returns>
15.         public static double NormSDist(double z)
16.         {
17.             return DistributionFunction(z);
18.         }
19.  
20.         /// <summary>
21.         /// The Probabilities density function.
22.         /// </summary>
23.         /// <param name="x">The x.</param>
24.         /// <returns></returns>
25.         public static double ProbabilityDensityFunction(double x)
26.         {
27.             const double ONE_OVER_SQRT_2_PI = 0.398942280401433;
28.             return ONE_OVER_SQRT_2_PI*Math.Exp(-0.5*x*x);
29.         }
30.  
31.         /// <summary>
32.         /// Returns the probability that a sample taken from the distribution lies inside the specified interval.
33.         /// </summary>
34.         /// <param name="minValue">The min value.</param>
35.         /// <param name="maxValue">The max value.</param>
36.         /// <returns>The probability that a sample taken from the distribution lies between minValue and maxValue</returns>
37.         public static double Probability(double minValue, double maxValue)
38.         {
39.             if (maxValue <= minValue) return 0.0;
40.             return DistributionFunction(maxValue) - DistributionFunction(minValue);
41.         }
42.  
43.         /// <summary>
44.         /// Calcula a distribuicao normal acumulada
45.         /// </summary>
46.         /// <param name="z">o fator de confiança z</param>
47.         /// <returns>a porcentagem acumulada até o valor de z [0; 1]</returns>        
48.         /// /// <remarks>
49.         /// Esta implementação segue HART 1968, conforme discutido em http://www4.ncsu.edu/unity/users/p/pfackler/www/ECG790C/accuratecumnorm.pdf
50.         /// que é o mesmo algoritmo que o Excel 2003, mas não o algoritmo do Excel XP, 2000 e anteriores.
51.         /// A precisão desta função é double, ao contrário da implementação anterior, que era single.
52.         /// </remarks>
53.         public static double DistributionFunction(double z)
54.         {
55.             if (double.IsNegativeInfinity(z))
56.             {
57.                 return 0.0;
58.             }
59.  
60.             if (double.IsPositiveInfinity(z))
61.             {
62.                 return 1.0;
63.             }
64.  
65.             if (z == 0.0)
66.             {
67.                 return 0.5;
68.             }
69.  
70.             double Cumnorm;
71.             double XAbs = Math.Abs(z);
72.             if (XAbs > 37.0)
73.             {
74.                 Cumnorm = 0.0;
75.             }
76.             else
77.             {
78.                 double Exponential = Math.Exp(-Math.Pow(XAbs, 2.0)/2.0);
79.                 if (XAbs < 7.07106781186547)
80.                 {
81.                     double build = 3.52624965998911E-02*XAbs + 0.700383064443688;
82.                     build = build*XAbs + 6.37396220353165;
83.                     build = build*XAbs + 33.912866078383;
84.                     build = build*XAbs + 112.079291497871;
85.                     build = build*XAbs + 221.213596169931;
86.                     build = build*XAbs + 220.206867912376;
87.                     Cumnorm = Exponential*build;
88.                     build = 8.83883476483184E-02*XAbs + 1.75566716318264;
89.                     build = build*XAbs + 16.064177579207;
90.                     build = build*XAbs + 86.7807322029461;
91.                     build = build*XAbs + 296.564248779674;
92.                     build = build*XAbs + 637.333633378831;
93.                     build = build*XAbs + 793.826512519948;
94.                     build = build*XAbs + 440.413735824752;
95.                     Cumnorm = Cumnorm/build;
96.                 }
97.                 else
98.                 {
99.                     double build = XAbs + 0.65;
100.                     build = XAbs + 4.0/build;
101.                     build = XAbs + 3.0/build;
102.                     build = XAbs + 2.0/build;
103.                     build = XAbs + 1.0/build;
104.                     Cumnorm = Exponential/build/2.506628274631;
105.                 }
106.             }
107.             if (z > 0.0)
108.             {
109.                 Cumnorm = 1.0 - Cumnorm;
110.             }
111.  
112.             return Cumnorm;
113.         }
114.  
115.         /// <summary>
116.         /// Retorna a função de sobrevivencia
117.         /// </summary>
118.         /// <param name="z">The z.</param>
119.         /// <returns></returns>
120.         public static double SurvivorDistributionFunction(double z)
121.         {
122.             return 1.0 - DistributionFunction(z);
123.         }
124.  
125.         /// <summary>
126.         /// Alias para a função StandardCumulativeInverse, com o mesmo nome que a função equivalente do Excel em inglês
127.         /// </summary>
128.         /// <param name="probability">probabilidade</param>
129.         /// <returns>O inverso da cumulativa padrão</returns>
130.         public static double NormSInv(double probability)
131.         {
132.             return InverseDistributionFunction(probability);
133.         }
134.  
135.         /// <summary>
136.         /// Calcula a inversa da cumulativa normal
137.         /// </summary>
138.         /// <param name="probability">probabilidade</param>
139.         /// <returns>O inverso da cumulativa padrão, 16 digitos de precisão</returns>
140.         /// <remarks>
141.         /// Retirado e adaptado do Fortran em 2010-10-26 a partir de http://lib.stat.cmu.edu/apstat/241
142.         ///
143.         ///    ALGORITHM AS241  APPL. STATIST. (1988) VOL. 37, NO. 3
144.         ///    Produces the normal deviate Z corresponding to a given lower
145.         ///    tail area of P; Z is accurate to about 1 part in 10**16.
146.         /// </remarks>
147.         public static double InverseDistributionFunction(double probability)
148.         {
149.             if ((probability < 0.0) || (probability > 1.0))
150.             {
151.                 throw new ArgumentOutOfRangeException("probability", probability, "Should be in [0; 1] interval.");
152.             }
153.  
154.             if (probability == 0.0)
155.             {
156.                 return double.NegativeInfinity;
157.             }
158.  
159.             if (probability == 1.0)
160.             {
161.                 return double.PositiveInfinity;
162.             }
163.  
164.             if (probability == 0.5)
165.             {
166.                 return 0.0;
167.             }
168.  
169.             const double SPLIT1 = 0.4250,
170.                          SPLIT2 = 5.0,
171.                          CONST1 = 0.1806250,
172.                          CONST2 = 1.60;
173.  
174.             //    Coefficients for P close to 0.5
175.             const double A0 = 3.38713287279636660800,
176.                          A1 = 1.3314166789178437745E+2,
177.                          A2 = 1.9715909503065514427E+3,
178.                          A3 = 1.3731693765509461125E+4,
179.                          A4 = 4.5921953931549871457E+4,
180.                          A5 = 6.7265770927008700853E+4,
181.                          A6 = 3.3430575583588128105E+4,
182.                          A7 = 2.5090809287301226727E+3,
183.                          B1 = 4.2313330701600911252E+1,
184.                          B2 = 6.8718700749205790830E+2,
185.                          B3 = 5.3941960214247511077E+3,
186.                          B4 = 2.1213794301586595867E+4,
187.                          B5 = 3.9307895800092710610E+4,
188.                          B6 = 2.8729085735721942674E+4,
189.                          B7 = 5.2264952788528545610E+3;
190.  
191.             //    Coefficients for P not close to 0, 0.5 or 1.
192.             const double C0 = 1.423437110749683577340,
193.                          C1 = 4.630337846156545295900,
194.                          C2 = 5.769497221460691405500,
195.                          C3 = 3.647848324763204605040,
196.                          C4 = 1.270458252452368382580,
197.                          C5 = 2.41780725177450611770E-1,
198.                          C6 = 2.27238449892691845833E-2,
199.                          C7 = 7.74545014278341407640E-4,
200.                          D1 = 2.053191626637758821870,
201.                          D2 = 1.676384830183803849400,
202.                          D3 = 6.89767334985100004550E-1,
203.                          D4 = 1.48103976427480074590E-1,
204.                          D5 = 1.51986665636164571966E-2,
205.                          D6 = 5.47593808499534494600E-4,
206.                          D7 = 1.05075007164441684324E-9;
207.  
208.             //    Coefficients for P near 0 or 1.
209.             const double E0 = 6.657904643501103777200,
210.                          E1 = 5.463784911164114369900,
211.                          E2 = 1.784826539917291335800,
212.                          E3 = 2.96560571828504891230E-1,
213.                          E4 = 2.65321895265761230930E-2,
214.                          E5 = 1.24266094738807843860E-3,
215.                          E6 = 2.71155556874348757815E-5,
216.                          E7 = 2.01033439929228813265E-7,
217.                          F1 = 5.99832206555887937690E-1,
218.                          F2 = 1.36929880922735805310E-1,
219.                          F3 = 1.48753612908506148525E-2,
220.                          F4 = 7.86869131145613259100E-4,
221.                          F5 = 1.84631831751005468180E-5,
222.                          F6 = 1.42151175831644588870E-7,
223.                          F7 = 2.04426310338993978564E-15;
224.  
225.             double r;
226.             double z;
227.             double q = probability - 0.5;
228.  
229.             if (Math.Abs(q) <= SPLIT1)
230.             {
231.                 r = CONST1 - q*q;
232.                 z = q*(((((((A7*r + A6)*r + A5)*r + A4)*r + A3)
233.                          *r + A2)*r + A1)*r + A0)/
234.                     (((((((B7*r + B6)*r + B5)*r + B4)*r + B3)
235.                        *r + B2)*r + B1)*r + 1.0);
236.                 return z;
237.             }
238.  
239.             if (q < 0.0)
240.             {
241.                 r = probability;
242.             }
243.             else
244.             {
245.                 r = 1.0 - probability;
246.             }
247.  
248.             if (r <= 0.0)
249.             {
250.                 throw new ArgumentOutOfRangeException("probability", probability, "Should be in [0; 1] interval.");
251.             }
252.  
253.             r = Math.Sqrt(-Math.Log(r));
254.             if (r <= SPLIT2)
255.             {
256.                 r = r - CONST2;
257.                 z = (((((((C7*r + C6)*r + C5)*r + C4)*r + C3)
258.                        *r + C2)*r + C1)*r + C0)/
259.                     (((((((D7*r + D6)*r + D5)*r + D4)*r + D3)
260.                        *r + D2)*r + D1)*r + 1.0);
261.             }
262.             else
263.             {
264.                 r = r - SPLIT2;
265.                 z = (((((((E7*r + E6)*r + E5)*r + E4)*r + E3)
266.                        *r + E2)*r + E1)*r + E0)/
267.                     (((((((F7*r + F6)*r + F5)*r + F4)*r + F3)
268.                        *r + F2)*r + F1)*r + 1.0);
269.             }
270.  
271.             if (q < 0.0)
272.             {
273.                 z = - z;
274.             }
275.  
276.             return z;
277.         }
278.     }

Os valores obtidos com NormSDist se comparam bem com os valores do Excel 2007, que parece usar algo parecido com Hart 1968, mas com outra escolha do ponto de quebra entre as funções. A diferença absoluta, no meio da distribuição (|z| < 3) é menor que 10^-15.  Nas bordas, com |z| >5 a diferença é menor que 10^-15, e por volta do ponto de quebra em |z| ~ 4 a diferença é menor que 4x10^-15.

Um teste interessante da implementação é computar p’ = NormSDist(NormSInv(p)) (com a classe acima) e comparar a diferença relativa entre p e p’ , d = (p’ – p) / p. Para p entre 0.01% e 99.99% o erro relativo absoluto máximo foi de 2.1x10^-14, conforme gráfico abaixo

Outro teste é computar z’ = NormSInv(NormSDist(z)) de modo similar.  Note que o erro aumenta bastante para valores de z extremos. De qualquer modo, se você estiver computando numa normal em lugares tão distantes do centro, é provável que esteja utilizando a distribuição incorreta, conforme West.

 

 

O uso da classe e dos algoritmos é livre, mas dê crédito a quem de crédito (Não, não para mim, só traduzi para C#). As consequências (especialmente as ruins) pelo uso, correto ou não, são suas, claro, mas se achar um Bug ou achar interessante uma mudança vou adorar saber.