您的问题完全不合适,正如其他人所指出的那样,所以我建议您先阅读一些内容并思考您真正想做什么。尽管如此,我可以提供一些指导,这实际上更多地涉及概念的梳理而非编程。基于问题和评论,您的问题应该重新表述为:“如何估计通过最大似然估计获得参数的AR(n)模型的参数协方差矩阵?”在这种情况下,我们可以回答如下。
- 为了讨论起见,让我们首先生成一些测试数据。
>> rng(0);
>> n = 10000;
>> x = rand(n,2);
2. 先通过运行OLS巧妙地获得初始条件。这是大多数人估计自回归的方式(向量自回归文献中的几乎全部),在某些条件下,OLS和MLE在渐近意义下是等价的。如果要包括常数项(如果您使用的是非去均值数据,则需要这样做),请包括一个零列。在输出中,b是您的参数估计向量,C是相应的标准误向量。
>> X = [x,ones(n,1)];
>> [b,C]=lscov(X,y)
b =
4.9825
9.9501
20.0227
C =
0.0347
0.0345
0.0266
在进行极大似然估计之前,您需要首先有一个似然函数。我们可以将其构建为几个简单的匿名函数的组合。我们还需要构建一个初始条件向量,其中包括预测误差的标准偏差,因为这是我们将使用MLE估计的内容之一。在这种情况下,我假设您希望误差服从正态分布,并基于此进行示例,但您没有说明,当然可以选择任何分布(不使用正态分布通常是不进行OLS的动机)。此外,我将构建一个负对数似然,因为我们正在使用最小化程序。
>> err = @(b) y - X*b
err =
@(b)y-X*b
>> std(err(b))
ans =
0.9998
>> b0 = [b; std(err(b))];
>> nll = @(b) (n/2)*log(2*pi*b(end)^2) + (1/(2*b(end)^2))*sum(err(b(1:end-1)).^2);
- 接下来使用某种最小化程序(例如,
fminsearch
,fminunc
等)进行MLE。
>> bmle = fminsearch(nll,b0)
bmle =
4.9825
9.9501
20.0227
0.9997
毫不意外的是,估计值几乎与我们在OLS下获得的相同,这正是我们在误差服从正态分布时所期望的,也是为什么大多数人选择只进行OLS,除非有充分的理由认为误差是非正态的。协方差矩阵可以通过得分的外积来估计,在正态情况下这是一个特别简单的表达式。
>> inv(X'*X/bmle(end))
ans =
0.0012 0.0000 -0.0006
0.0000 0.0012 -0.0006
-0.0006 -0.0006 0.0007
最后,标准误差与我们在最小二乘法中得到的结果相匹配。
>> sqrt(diag(inv(X'*X/bmle(end))))
ans =
0.0347
0.0345
0.0266
编辑:抱歉,我刚意识到我的测试数据是横向截面数据而不是时间序列数据。我会尽快修复它。但是用于估计模型的方法保持不变。
cov
的东西,但不知道AR。你有什么AR模型的输出?cov()
函数是用于在给定变量的一堆观测/样本(输入=观测x变量数组)的情况下数值计算协方差的,不确定这里是否适用。如果您想要系数的协方差,是否有一个公式来产生系数本身的样本?(a_i
值是否为这些系数?) AR可能有一个封闭形式的解,但我不知道它是什么。 - Andrew Jankecov
函数并不是你要找的。你应该使用估计的系数a_j
构建p
-by-p
矩阵A_1
和A_2
,因为ML估计器还提供了噪声方差sigma
的估计值(假设它确实提供)。 - lmillefiori