标签:
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()
标签:
原文地址:http://www.cnblogs.com/prettysmc/p/5500070.html