[游戏开发]2021-11-08

TAR模型（门限自回归模型）

TAR模型的实质是分段的AR模型，，它的基本思路是，在观测时序{$x_i$}的取值范围内引入$l?1$个门限值$r_i$($i$=1,2,…$l?1$),将该范围分成$l$个区间，可用$r_0$,$r_l$分别表示上界和下界，并根据延迟步数$d$将{$x_i$}按{$x_{i?d}$}值的大小分配到不同的门限区间内，再对区间内的{$x_i$}采用不同的自回归模型（AR模型），从而形成时间序列的非线性动态描述，其模型形式如下：
$x_t$=${\Phi }_0^j$+${\sum }_{i=1}^{P_j}$${\Phi }_i^j$$X_{t?i}$+${\epsilon }_t^j$,$r_j$,$j$=1,2,…$l$

式中${\epsilon }_t^j$是第$l$个相互独立的正态分布白噪声序列，$d$为延迟步数（非负整数），$r_j$为门限值，$l$为门限区间的个数，${\Phi }_0^j$为第$j$个门限区间的自回归系数，$P_j$为第$j$个门限区间AR模型的阶数。
由于TAR模型实质是区分区间的AR模型，建模时刻沿用AR模型的参数估计方法模型检验准则，其建模实质是第一个对$d,l,r_j,p_j$,${\Phi }_i^j$的多维寻优问题，可在给定哥哥参数取值范围的条件下，使用遗传算法、网格搜索等方法搜索最有参数。门限区间个数$l$的选取，理论上可以选取若干个，但在实际应用中往往选择1对就可以满足，因此通常将门限区间个数$l$取为2。假设用于训练的数据样本量为$n$,则可依次取0.3$n$，0.4$n$，0.5$n$，0.6$n$，0.7$n$，的对应的值作为门限候选值。对于候选值中的任意门限值$r_1$，设定最大门限延迟量为$d_{max}$，对于任意的$1<=d<=d_{max}$，我们可以将时序数据分成两类，一类时序值其$d$阶延迟值小于等于$r_1$，另一类时序值大于$r_1$。针对这爱玲组数据，分别建立AR模型，并计算出这艾灵格模型对应的AIC值，当他们的和最小时，对应的$d,r_1,p_1,p_2$,${\Phi }_j^1$($i=0,2,...p_1$),${\Phi }_i^2$($i=0,2,...,p_2$)参数即为所求，其中AIC的公式计算如下：
$AI{C}_j$ = $n_{j}ln$(${\sigma }_j^2$)+2($p_j$+2),$j=1,2$
式中$n_j$为第$j$个门限区间的样本个数，${\sigma }_j^2$为第$j$个门限区间的阿姨你根本残差的方差，$p_j$为对应的自回归模型阶数。

TAR适用于周期性波动、非线性影响等情况。

R语言
vdata = read.table("f:\\sunspot-1700-2010.txt",header=F)
dim(vdata)
#使用未来50年的数据进行验证，以前的数据用于建模
modelData = vdata$V2[1:261]
n <- length(modelData)
#设置候选门限值
thresholdV = sort(modelData)[round((seq(from=30,to=70,by=3)/100)*n)]
#设置门限延迟量dmax、自回归最大阶数，默认最小AIC值
dmax = 5
pmax = 5
minAIC = le+10
#在指定们线延迟量、结束及门限值的前提下，返回对应自回归AIC值和自回归系数
getModelInfo <- function(tsoobj,d,p,r,isUp=True){
if(isUp){
dstSet = which(tsobj>r)+d
}else{
dstSet = which(tsobj<=r)+d}
tmpdata = NULL
#重建基础数据集
#xt = a0+a1*x(t-1)+...+ap*x(t-p)
for(i in dstSet){
if (i>p &&i<length(tsobj)){
tmpdata = rbind(tmpdata,tspbj[i:(i-p)])
}
}
x = cbind(1,tmpdata[,2:(p+1)])
coef = solve(t(x)%*%x)%*%tmpdata[,1]
epsion = tmpdata[,1]-x%*%coef
aic = drop(nrow(tmpdata)*log(var(epsion))+2*(p+2))
return (list(aic = aic,coef = coef))
}
#选择最优参数
for(tsv in thresholdV){
for (d in 1:dmax){
for (p1 in 1:pmax){#<=r1
model1 = getModelInfo(modelData,d ,p = p1,r = tsv,isUp = F)
for (p2 in 1:pmax)#>r1
{model2 = getModelInfo(modelData,d ,p = p1,r = tsv,isUp = T)
if ((modell$aic+model2$aic)<minAIC)
{
minAIC = (model1$aic+model2$aic)
a_tsv = tsv
a_d = d
a_p1 = p1
a_p2 = p2
coef1 = model1$coef
coef2 = model2$coef
print(minAIC)
}
}
}
}
}

#使用求出的参数，对未来50年的数据逐年预测
predsData = NULL
for(i in 262:311)
{
t0 = vdata$V2[(i-a_d)]
if(t0 <= a_tsv)
{
predsData=c(predsData,sum(c(1,vdata$V2[(i-1):(i-a_p1)])*coef1))
}else{
predsData = c(predsData,sum(c(1,vdata$V2[(i-1):(i-a_p2)])*coef2)}
}


#调包
library（TSA）
vdata = read.table("f:\\sunspot-1700-2010.txt",header=F)
#使用未来50年的数据进行验证，以前的数据用于建模
篇reds DATa= NULL
for(i in 1:50)
{
modelData = vdata$V2[1:(261+i-1)]
tar.obj = tar(y=modelData,p1=5,p2=3,d=3,a=3,b=7)
predsData = c(predsData,predict(tar.obj,n.ahead=1)$fit)
}