Sugestão de metodologia para cálculo de Intervalos de Confiança

Conforme mencionado na LDO de 2018, o modelo oficial do governo se define como determinístico:

“[...] ou seja, a partir da fixação de um conjunto de variáveis, o modelo determina de maneira única seus resultados [...]

Como se trabalha com probabilidades, não necessariamente todos os eventos previstos podem acontecer. O modelo da LDO é determinístico por trabalhar apenas com médias (ex: média de pessoas que se aposentarão) e não considera diferentes cenários onde isso pode não ocorrer, ou seja, situações diferentes do comportamento médio.

Este documento busca apresentar uma forma diferente de se projetar estoques considerando diferentes cenários onde nem sempre os segurados irão se aposentar.


In [ ]:
import sys
sys.path.insert(0, '../')

from util.tabelas import LerTabelas
import modelos.fazenda as fz
%matplotlib inline


###################### Parâmetros de simulação ###############################

# Período de projeção 
periodo = list(range(2015, 2061))

# Ano de referência para cálculo das probabilidades
ano_probabilidade = 2014

#############################################################################

arquivo = '../dados/dados_fazenda.xlsx'

# Abri o arquivo
dados = LerTabelas(arquivo)

# Beneficio que será avaliado
id_aptcn = 'AtcnUrbPisoH'

# Lista de Ids dos beneficios
ids_estoques = dados.get_id_beneficios([], 'Es')
ids_concessoes = dados.get_id_beneficios([], 'Co')
ids_cessacoes = dados.get_id_beneficios([], 'Ce')

# Obtem as tabelas e armazena nos dicionários correspondentes
estoques = dados.get_tabelas(ids_estoques)
concessoes = dados.get_tabelas(ids_concessoes)
cessacoes = dados.get_tabelas(ids_cessacoes)
populacao = dados.get_tabelas(dados.ids_pop_ibge)
populacao_pnad = dados.get_tabelas(dados.ids_pop_pnad)

In [ ]:
# Calcula taxas de urbanização, participação e ocupação
taxas = fz.calc_taxas(populacao_pnad, periodo)

# Calcula: Pop Urbana|Rural, PEA e Pop Ocupada, Contribuintes e Segurados Rurais
segurados = fz.calc_demografia(populacao, taxas)

# Calcula as probabilidades de entrada em benefício e morte
probabilidades = fz.calc_probabilidades(populacao, segurados, estoques,
                                     concessoes, cessacoes, periodo)

In [ ]:
probabilidades[id_aptcn][2014].plot();

In [ ]:
import numpy as np

n_execucoes = 15

def calc_qtd_apos(num_segurados, prob_conc, seed):
        
    # Define o seed para geração de números aleatórios
    np.random.seed(seed)
            
    # Gera números aletórios entre 0 e 1 de acordo com a quantidade de segurados
    chance_segurados = np.random.rand(num_segurados)
    
    concessoes = 0

    # Determina quantos irão receber o benefício 
    for chance in chance_segurados:
        # calcula a probabilidade de cada um dos segurados receber o benefício
        # Se o número aleatório for menor ou igual a probabilidade o segurado irá se aposentar
        if chance <= prob_conc: 
            concessoes += 1
            
    #print('{} - {} - {} - {}'.format(num_segurados, prob_conc, chance_segurados[0], concessoes))
    #print(chance_segurados)
    
    return concessoes


# Cria a função que calcula o estoque usando as equações da LDO
def calc_estoq_apos(est, conc, prob, seg, periodo):
    
    # Calcula somente para a Aposentadoria por Tempo de Contribuição
    ids_apos= ['AtcnUrbPisoH']

    # Cria o objeto dados que possui os IDs das tabelas
    dados = LerTabelas()

    # Obtem as aposentadorias para todas as clientelas e sexos
    lista_benef = dados.get_id_beneficios(ids_apos)
        
    for benef in lista_benef:
        # Verifica se o beneficio existe no Estoque
        if benef in est:
        
            sexo = benef[-1]                # Obtém o Sexo
            id_prob_morte = 'Mort'+ sexo    # ex: MortH
            id_fam = 'fam'+benef            # fator de ajuste de mortalidade            
            id_segurado = dados.get_id_segurados(benef)  # ex: CsmUrbH
            
            for ano in periodo:                
                # Adiciona uma nova coluna (ano) no DataFrame com valores zero
                est[benef][ano] = 0
                
                # 1 a 90 anos - Equação 11 da LDO de 2018
                for idade in range(1,91): 
                    est_ano_anterior = est[benef][ano-1][idade-1]
                    prob_sobreviver = 1 - prob[id_prob_morte][ano][idade] * prob[id_fam][ano][idade]
                    entradas = seg[id_segurado][ano][idade] * prob[benef][ano][idade]
                    
                    # Eq. 11
                    est[benef].loc[idade, ano] = est_ano_anterior * prob_sobreviver + entradas     # Eq. 11
                    # Salva a quantidade de concessões para uso posterior
                    conc[benef].loc[idade,ano] = entradas
                
                # Calculo para a idade zero
                est[benef].loc[0, ano] = seg[id_segurado][ano][0] * prob[benef][ano][0]
                # Salva a quantidade de concessões para uso posterior
                conc[benef].loc[0, ano] = est[benef].loc[0, ano]
                
                # Ajuste para a idade de 90+ anos (modelo UFPA) - REVISAR
                est[benef].loc[90, ano] = est[benef].loc[90, ano] + est[benef].loc[90, ano - 1]
                

    return est



# Cria a função que calcula o estoque usando as equações da LDO
def calc_estoq_apos_estocastico(est_2, conc_2, prob, seg, periodo):
    
    # Calcula somente para a Aposentadoria por Tempo de Contribuição
    ids_apos= ['AtcnUrbPisoH']

    # Cria o objeto dados que possui os IDs das tabelas
    dados = LerTabelas()

    # Obtem as aposentadorias para todas as clientelas e sexos
    lista_benef = dados.get_id_beneficios(ids_apos)
        
    for benef in lista_benef:
        # Verifica se o beneficio existe no Estoque
        if benef in est_2:
        
            sexo = benef[-1]                # Obtém o Sexo
            id_prob_morte = 'Mort'+ sexo    # ex: MortH
            id_fam = 'fam'+benef            # fator de ajuste de mortalidade            
            id_segurado = dados.get_id_segurados(benef)  # ex: CsmUrbH
            
            est_seed = dict()
            conc_seed = dict()
            
            seeds = range(0,n_execucoes)
            
            # Para cada um dos seeds
            for seed in seeds:
            
                estoq = copy.deepcopy(est_2)
                conces = copy.deepcopy(conc_2)
                                
                for ano in periodo:                
                    # Adiciona uma nova coluna (ano) no DataFrame com valores zero
                    estoq[benef][ano] = 0

                    # 1 a 90 anos - Equação 11 da LDO de 2018
                    for idade in range(1,91): 
                        est_ano_anterior = estoq[benef][ano-1][idade-1]
                        prob_sobreviver = 1 - prob[id_prob_morte][ano][idade] * prob[id_fam][ano][idade]

                        # MUDEI AQUI
                        qtd_segurados = seg[id_segurado].loc[idade, ano]
                        prob_conc = prob[benef].loc[idade, ano]

                        #print('{} - {}'.format(ano, idade))

                        if qtd_segurados < 1 or prob_conc == 0.0:
                            entradas = 0
                        else:                        
                            entradas = calc_qtd_apos(qtd_segurados, prob_conc, seed)                        
                            #print('{} - {} - {} - D = {} - E = {}'.format(ano, idade, qtd_segurados, entradas))

                        # Eq. 11
                        estoq[benef].loc[idade, ano] = est_ano_anterior * prob_sobreviver + entradas     # Eq. 11
                        # Salva a quantidade de concessões para uso posterior
                        conces[benef].loc[idade,ano] = entradas

                    # Calculo para a idade zero
                    estoq[benef].loc[0, ano] = seg[id_segurado][ano][0] * prob[benef][ano][0]
                    # Salva a quantidade de concessões para uso posterior
                    conces[benef].loc[0, ano] = estoq[benef].loc[0, ano]

                    # Ajuste para a idade de 90+ anos (modelo UFPA) - REVISAR
                    estoq[benef].loc[90, ano] = estoq[benef].loc[90, ano] + estoq[benef].loc[90, ano - 1]
                
                est_seed[seed] = estoq
                conc_seed[seed] = conces
                
        return est_seed, conc_seed


# Cria uma cópia dos dados de estoque
import copy
estoques_2 = copy.deepcopy(estoques)
concessoes_2 = copy.deepcopy(concessoes)

# Executa a função
est_deterministico = calc_estoq_apos(estoques, concessoes, probabilidades, segurados, periodo)
est_estocastico, conc_estocastico = calc_estoq_apos_estocastico(estoques_2, concessoes_2, probabilidades, segurados, periodo)

In [ ]:
concessoes[id_aptcn][2050].plot();
for seed in conc_estocastico.keys():
    conc_estocastico[seed][id_aptcn][2050].plot();

In [ ]:
import matplotlib.pyplot as plt

fig, ax = plt.subplots()
ax.plot(concessoes[id_aptcn].sum().loc[2018:], '-', linewidth=2, label='Media')

for seed in conc_estocastico.keys():
    ax.plot(conc_estocastico[seed][id_aptcn].sum().loc[2018:], ':', linewidth=2, label=seed)    
ax.set_ylabel('Número de Concessões')
ax.set_xlabel('Anos')
ax.set_title('Cálculo das Concessões para diferentes seeds')
ax.legend(loc='center left', bbox_to_anchor=(1, 0.5))
plt.show()

In [ ]:
import numpy as np
import scipy as sp
import scipy.stats
import pandas as pd


def mean_confidence_interval(data, confidence=0.95):
    a = 1.0*np.array(data)
    n = len(a)
    m, se = np.mean(a), scipy.stats.sem(a)
    h = se * sp.stats.t._ppf((1+confidence)/2., n-1)
    return m, m-h, m+h

# Calculando intervalo de confiança das concessoes
resultado = pd.DataFrame(index=periodo, columns=['media', 'LimSup', 'LimInf'])

for ano in periodo:
    temp = []  # lista para salvar os valores de concessões
    for seed in conc_estocastico.keys():
        # Sala o total de concessoes na lista
        temp.append(conc_estocastico[seed][id_aptcn].sum().loc[ano])
        
    #print('{} - {}'.format(ano, temp))
    resultado['media'][ano], resultado['LimSup'][ano], resultado['LimInf'][ano] = mean_confidence_interval(temp)

In [ ]:
#Testes

from math import sqrt

print(temp)

print('LimInf - Média - LimSup')

##### Metodo 01
med, hi, lw= mean_confidence_interval(temp)
print('{} - {} - {}'.format(hi, med, lw))

##### Método 02
t2 = 1.0*np.array(temp)
med2 = t2.mean()
std2 = t2.std(ddof=1)
d = 1.96 * (std2/sqrt(len(t2)))

print('{} - {} - {}'.format(med2-d, med2, med2+d))

##### Método 03
from scipy import stats

mean, sigma = np.mean(t2), np.std(t2)
conf_int = stats.norm.interval(0.95, loc=mean, scale=sigma)
print('{} - {} - {}'.format(conf_int[0], mean, conf_int[1]))

In [ ]:
fig, ax = plt.subplots()
ax.plot(resultado['media'], '-', linewidth=2)
#ax.plot(resultado['LimSup'], '-', linewidth=2, color=(0, 0, 0))
#ax.plot(resultado['LimInf'], '-', linewidth=2, color=(0, 0, 0))
c = 1
ax.fill_between(periodo, resultado['LimSup'].tolist(), resultado['LimInf'].tolist(),
                     color=(0.8-0.20*c, 0.8-0.20*c, 1.0),     
                     facecolor=(0.8-0.20*c, 0.8-0.20*c, 1.0),
                     label='95%')     

ax.set_ylabel('Concessões')
ax.set_xlabel('Anos')
ax.set_title('Projeção das Concessões para Aposentadoria por TC')
ax.legend(loc='center left', bbox_to_anchor=(1, 0.5))
plt.grid(True)
plt.show()

In [4]:
import time
import sys

for x in range(101):
    print("\rProgress {}%".format(x), end="")
    time.sleep(0.1)
    sys.stdout.flush()


Progress 100%

In [79]:
import numpy as np

def calc_qtd_apos(num_segurados, prob_conc, seed):
        
    # Define o seed para geração de números aleatórios
    np.random.seed(seed)
            
    # Gera números aletórios entre 0 e 1 de acordo com a quantidade de segurados
    chance_segurados = np.random.rand(num_segurados)
    
    concessoes = 0

    # Determina quantos irão receber o benefício 
    for chance in chance_segurados:
        # calcula a probabilidade de cada um dos segurados receber o benefício
        # Se o número aleatório for menor ou igual a probabilidade o segurado irá se aposentar
        if chance <= prob_conc: 
            concessoes += 1
            
    #print('{} - {} - {} - {}'.format(num_segurados, prob_conc, chance_segurados[0], concessoes))
    #print(chance_segurados)
    
    return concessoes

print(calc_qtd_apos(1000, 0.25, 1))
print(calc_qtd_apos(1000, 0.25, 2))
print(calc_qtd_apos(1000, 0.25, 3))
print(calc_qtd_apos(1000, 0.25, 4))
print(calc_qtd_apos(1000, 0.25, 5))
print()

np.random.seed(1)
#for i in range(30):
print(np.mean(np.random.binomial(1000, 0.25,10000)))


249
256
236
235
247

249.7798