INAR模型参数估计的计算机实现

Posted 2020-07-06

1. 估计方法：Yule-Walker估计

.

2. 运行环境spyder2(python科学计算编辑器)

3. 实现代码

# -*- coding: utf-8 -*-
import numpy as np#创建数组用的包
#设置初值
A=9

f = open(‘/Users/gaojiaxing/desktop/data.txt‘,‘a‘,encoding = ‘utf-8‘)
f.write(‘时间跨度‘+‘ ‘+‘p值‘+‘ ‘+‘lamb值‘+‘ ‘+‘p均方误差‘+‘ ‘+‘lamb均方误差‘+‘ ‘+‘p的估计‘+‘ ‘+‘p的估计‘+‘ ‘+‘p的估计‘+‘ ‘+‘lamb的估计‘+‘ ‘+‘lamb的估计‘+‘ ‘+‘lamb的估计‘+‘\\n‘)
f.close()

def sample_cov(alist,k):
    ‘‘‘
    这个函数是用来计算样本的自协方差函数
    ‘‘‘
    n = len(alist)
    sum1 = 0
    for i in range(n-k):
        sum1 = sum1 + (alist[i]-np.mean(alist))*(alist[i+k]-np.mean(alist))
    return sum1/n
    
def calc_lamb(alist):
    ‘‘‘
    yule-walker估计求误差项的
    ‘‘‘
    n = len(alist)
    error_sum = 0
    for j in range(n-1):
        error_sum = error_sum + alist[j+1]-p*alist[j]
    return error_sum
  
for time_len in range(50,250,50):#时间序列长度
    for p in np.arange(0,1,0.1):#让p值取0,0.1，0.2，...0.9的所有结果
        for lamb in np.arange(0.2,1.4,0.4):#lamb取值0.2，0.6，1
            
            claim_time_series = []#模拟的最终1000组时间序列的存放在这里
    
            for i in range(1000):
                one_case = list(range(time_len))#one_case是一个理赔时间序列
                one_case[0] = int(lamb/(1-p))
                for j in range(1,time_len):
                    one_case[j] = (np.random.binomial(one_case[j-1],p,size=1)+                    np.random.poisson(lamb*np.random.gamma(shape=A,scale=1/A,size=1)))[0]
                    
                claim_time_series.append(one_case)
                
            p_hat = []#根据样本估计出来的P估计值将要存放的地方
            for j in range(1000):
                p_hat.append(sample_cov(claim_time_series[j],1)/sample_cov(claim_time_series[j],0))
                
            p_array = np.array(p_hat).reshape(1000)#转换成数组，便于运算

            p_mse = sum((p_array-p*np.ones(1000))**2)/1000#均方误差

            lamb_hat = []
            for i in range(1000):
                lamb_hat.append(calc_lamb(claim_time_series[i])/(time_len-1))
                
            lamb_array = np.array(lamb_hat).reshape(1000)
            
            lamb_mse = sum((lamb_array-lamb*np.ones(1000))**2)/1000

            f = open(‘/Users/gaojiaxing/desktop/data.txt‘,‘a‘,encoding = ‘utf-8‘)
            f.write(str(time_len)+‘ ‘+str(p)+‘ ‘+str(lamb)+‘ ‘+str(p_mse)+‘ ‘+            str(lamb_mse)+‘ ‘+str(p_array[0])+‘ ‘+str(p_array[1])+‘ ‘+str(p_array[2])+‘ ‘+            str(lamb_array[0])+‘ ‘+str(lamb_array[1])+‘ ‘+str(lamb_array[2])+‘\\n‘)
            f.close()