1、时间序列分析产生序列的matlab代码见最后一、AR(p)序列模型（一）产生数据先用matlab产生AR（4）模型的数据 X(t)=a(1)X(t-1)+ a(2)X(t-2)+ a(3)X(t-3)+ a(4)X(t-4)+(t) , (t)服从N(0,2).先选取四个多项式单位根，z1,z2,z3,z4,他们的模都大于一，然后产生相应的多项式，得到一组系数a(1),a(2),a(3),a(4)，满足模型序列是平稳列。 p=4; r=2,3,4,5; %模型特征多项式的根下面产生该模型相应的数据,500个点；(t)服从N(0,c). N=500; c=1; X,a=generate_ar(p
2、,r,N,c);%a是模型系数模型为：X(t)= 1.2833*X(t-1)-0.5917*X(t-2)+0.1167*X(t-3)-0.0083*X(t-4)+(t) 产生的序列图如下： plot(X) xlabel(t); ylabel(X(t);（二）单位根检验通过检验，序列没有单位根，没有常数项，所以是零均值平稳列，与原模型相符。（三）给模型定阶，求系数从自相关偏相关图来看，序列对PACF截断，对ACF衰减，又因为序列是零均值平稳列，暂定为AR（3）。虽然从上面已知该序列是零均值平稳列，先加上常数项试试，做个检验。常数项系数不是很显著。AR（3）系数特别不显著，与原模型有差。先去掉AR
3、（3）试试。常数项仍旧不显著，这与原模型一致。去掉试试这个模型的R2与AIC跟别的差不多，各个系数均显著，较别的模型也较为简单，所以先记上面这个模型为模型1,即X(t)=1.250384X(t-1)-0.52247X(t-2)下面再加上AR（3）试试；仍旧不显著，鉴于已知原模型，所以加上AR（4）看看效果记上模型为模型2：X(t)=1.253088X(t-1)-0.5037X(t-2)- 0.059X(t-3)+0.0556X(t-4)AR(3),AR(4)仍旧不显著，进而在不知道原模型的情况下，最佳模型应该是模型1：X(t)=1.250384X(t-1)-0.52247X(t-2)拟合图如下
4、：看着效果已经很好了，下面截一小段看看： 虽然与原模型不完全一致，但是看拟合效果，已经很不错了，而且从原模型也可以看到，其实AR(3),AR(4)的系数本来就很小，所以建模时没有也正常，如果加上AR(3),AR(4)，也即模型2：X(t)=1.253088X(t-1)-0.5037X(t-2)- 0.059X(t-3)+0.0556X(t-4)原模型：X(t)= 1.2833X(t-1)-0.5917X(t-2)+0.1167X(t-3)+ 0.008X(t-4)+(t)两个模型在AR(3)上相差较大。下图为模型二的拟合图。上面两个图是模型一和模型二的拟合对比图，两个模型拟合效果差不多，几乎重
5、合，本来模型2中ar(3),ar(4)系数较小，所以影响也不大。看原模型与真值的差异，它是由随机项白噪声产生的：从上面的图可以看到加入随机项后，有时候新建的模型要比原模型的拟合的更好一些，还比原模型简单，也没有滞后效果。再与简单模型X(t)=X(t-1)（记为模型3）做对比感觉模型1要比模型3整体稍微好一点，不过差不太多，而且原模型的AR(1)系数为1.25，相邻的数据相关性较大，所以它们差别不大，用昨日预测今日也是比较合理的，再者，取值范围较小，所以误差也不会太大，差别也不是很大。下面是他们两残差图的比较：系列1是模型1的残差，系列2是模型3的残差，感觉模型1还是要比模型3即简单模型要好一些
6、，但也差不太多。所以最佳模型为：X(t)=1.250384X(t-1)-0.52247X(t-2)。二、ARMA（p,q）模型（一） 产生数据ARMA(5,3),即p=5,q=3.先随意产生一组使得模型序列是平稳序列的向量a1,a5；b1,b3；b0=1;因为要满足是平稳列，需要特征多项式A的根1，特征多项式B的根1；且它两之间没有相同的根，所以取A的根 r1=1.8 ,-1.1 ,-1.6 ,1.7, -1.4;B的根r2=1.5, 2.1 ,1.3; r1=1.8 ,-1.1 ,-1.6 ,1.7, -1.4; r2=1.5,2.1,1.3; p=5;q=3; c=20;%常数项 N=50
7、0; X,a,b=generate_arma(p,q,r1,r2,N,c);经过多次尝试，尽量选了a(5),b(3)系数较为显著的一个模型，因为上面AR(4)的时候，在建模的时候ar(4)一直不通过，但是原模型有ar(4),只是它的系数很小，这次正好检验一下是不是因为原模型系数小，所以不通过。此外，也尽量选择ar(1)是负的，这样的话用前一天来预测今天效果是不就不好了。而且这次给模型加上了常数项。所以产生的模型为：X(t)=20-1.1046X(t-1)+0.5809X(t-2)+0.7626X(t-3)-0.0796X(t-4)-0.1326X(t-5)-1.9121(t-1)+1.1966
8、(t-2)-0.2442(t-3)+(t); plot(X) title(arma(5,3) 这次这个序列与别比就奇葩多了，波动相当大，下面看建模效果如何。（一） 平稳性检验通过检验，序列为平稳列，有常数项，不是零均值的，与原模型一致。（二）、建模、定阶、求系数序列对PACF截断，对ACF呈指数，从这个来看，模型应该是AR序列，先顺着这个趋势看看。AIC=2.841，AR(9)，AR(10)系数稍有不显著，去掉后如下，AIC比前面稍小一些，这个可以近似认为各系数均显著，下图为拟合图残差在几乎都在0-1之间，拟合效果已经很好了，把这个模型记作模型1：X(t)=20.54676-2.9484X(t
9、-1)-3.72X(t-2)-3.04X(t-3)-2.27X(t-4)-1.74X(t-5)-1. 2X(t-6)-0.5735X(t-7)-0.1271X(t-8)下面去掉ar(8)看有什么效果。这时，本来随着项数的减小，R2应该也会变小才对，但是反而增大，而且AIC也突然增加好多，再看拟合图：从坐标系明显看出，这个残差大一些。最后看AR（1）模型AIC很大，残差也较前面大很多，在0-5之间，前面是0-1之间。显然效果不好。 如果考虑它是ARMA模型，则走下路线把系数不显著的均去掉此时各系数几乎都显著，AIC也很小，下面是拟合图，相比较而言，如果是ARMA模型，综合考虑，上面这个已经是最好
10、的了记为模型2：X(t)=20.54676-2.11032(t-1)-1.25885X(t-2)+0.1375X(t-4)-0.8328(t-1)因为已知原模型是ARMA(5,3),所以从这个出发再试试。这个模型也很不错，R2,AIC都很小，有一个系数稍有不显著。由于随机项，各个系数与原模型系数差的较远。所以这三个模型都很好，残差几乎在0-1之间，相比于原值，几乎看不出差别，不过模型2最简洁。function X,a=generate_ar(p,r,N,c) %a为模型系数，X产生的数据。a1=poly(r);I=1:p;a(I)=-a1(5-I)/a1(p+1);Y=zeros(4,1);e
11、ps=1:500;for t=5:N+55 eps(t)=c*randn(1); Y(t)=sum(a.*Y(t-I)+eps(t); %Y(t)=a(1)*Y(t-1)+a(2)*Y(t-2)+a(3)*Y(t-3)+a(4)*Y(t-4)+eps(t);endJ=1:500;X(J)=Y(J+55);function X,a,b=generate_arma(p,q,r1,r2,N,c) %a,b为模型系数，c是常数项=常数项，X产生的数据。a1=poly(r1);b1=poly(r2);I=1:p;Iq=1:q;a(I)=-a1(p+1-I)/a1(p+1);b(Iq)=b1(q+1-Iq
12、)/b1(q+1);Y=zeros(p,1);eps=1:N;eps(I)=randn(1,p);for t=p+1:N+p+1+300 eps(t)=randn(1); Y(t)=c+sum(a.*Y(t-I)+sum(b.*eps(t-Iq)+eps(t); %Y(t)=a(1)*Y(t-1)+a(2)*Y(t-2)+a(3)*Y(t-3)+a(4)*Y(t-4)+eps(t);endJ=1:N;X(J)=Y(J+p+301);function Y=moni(y,p,q,a,b,N) c=1;%a为模型系数，X产生的数据。I=1:p;Iq=1:q;Y(I)=y(I);eps=1:N;eps
13、(I)=c*randn(1,p);for t=p+1:N eps(t)=c*randn(1); Y(t)=sum(a.*y(t-I)+sum(b.*eps(t-Iq); %Y(t)=a(1)*Y(t-1)+a(2)*Y(t-2)+a(3)*Y(t-3)+a(4)*Y(t-4)+eps(t);endfunction Y=moni1(y,p,a1,N) c=1;%a为模型系数，X产生的数据。I=1:p;Y(I)=y(I);I2=p+1:N;for t=p+1:N Y(t)=sum(a1.*y(t-I); %Y(t)=a(1)*Y(t-1)+a(2)*Y(t-2)+a(3)*Y(t-3)+a(4)*
14、Y(t-4)+eps(t);end附录资料：MATLAB函数和命令的用法Binocdf二项式累积分布函数语法格式Y = binocdf(X,N,P)函数功能Y = binocdf(X,N,P) 计算X中每个X(i)的二项式累积分布函数，其中，N中对应的N(i)为试验数，P中对应的P(i)为每次试验成功的概率。Y, N, 和 P 的大小类型相同，可以是向量、矩阵或多维数组。输入的标量将扩展成一个数组，使其大小类型与其它输入相一致。The values in N must all be positive integers, the values in X must lie on the inter
15、val 0,N, and the values in P must lie on the interval 0, 1.The binomial cdf for a given value x and a given pair of parameters n and p is The result, y, is the probability of observing up to x successes in n independent trials, where the probability of success in any given trial is p. The indicator
16、function I(0,1,.,n)(i)ensures that x only adopts values of 0,1,.,n.示例若一个棒球队在一个赛季要比赛162场，每场比赛取胜的机会是50-50，则该队取胜超过100 场的概率为： 1-binocdf(100,162,0.5)ans = 0.0010433相关函数binofit | binoinv | binopdf | binornd | binostat | cdf附：二项式分布(binomial distribution )定义二项分布的概率密度函数为where k is the number of successes in
17、n trials of a Bernoulli process with probability of success p.The binomial distribution is discrete, defined for integers k = 0, 1, 2, . n, where it is nonzero.背景The binomial distribution models the total number of successes in repeated trials from an infinite population under the following conditio
18、ns:Only two outcomes are possible on each of ntrials.The probability of success for each trial is constant.All trials are independent of each other.The binomial distribution is a generalization of the Bernoulli distribution; it generalizes to the multinomial distribution.参数Suppose you are collecting
19、 data from a widget manufacturing process, and you record the number of widgets within specification in each batch of100. You might be interested in the probability that an individual widget is within specification. Parameter estimation is the process of determining the parameter, p, of the binomial
20、 distribution that fits this data best in some sense.One popular criterion of goodness is to maximize the likelihood function. The likelihood has the same form as the binomial pdf above. But for the pdf, the parameters (nandp) are known constants and the variable isx. The likelihood function reverse
21、s the roles of the variables. Here, the sample values (the xs) are already observed. So they are the fixed constants. The variables are the unknown parameters. MLE involves calculating the value of p that give the highest likelihood given the particular set of data.The function binofit returns the M
22、LEs and confidence intervals for the parameters of the binomial distribution. Here is an example using random numbers from the binomial distribution with n=100 and p=0.9. r = binornd(100,0.9)r = 85 phat, pci = binofit(r,100)phat = 0.85pci = 0.76469 0.91355The MLE for parameterp is0.8800, compared to
23、 the true value of0.9. The 95% confidence interval forp goes from 0.7998 to0.9364, which includes the true value. In this made-up example you know the true value ofp. In experimentation you do not.示例The following commands generate a plot of the binomial pdf for n = 10 and p = 1/2.x = 0:10;y = binopd
24、f(x,10,0.5);plot(x,y,+) 相关内容Discrete Distributions附：二项式分布（网上）定义若某事件概率为p，现重复试验n次，该事件发生k次的概率为:P=C(k,n)pk(1-p)(n-k)C(k,n)表示组合数，即从n个事物中拿出k个的方法数。二项分布的概念考虑只有两种可能结果的随机试验，当成功的概率（）是恒定的，且各次试验相互独立，这种试验在统计学上称为贝努里试验（Bernoullitrial）。如果进行n次贝努里试验，取得成功次数为X（X=0,1,n）的概率可用下面的二项分布概率公式来描述：P=C(X,n) X(1-)(n-X)式中的n为独立的贝努里试验
25、次数，为成功的概率，（1-）为失败的概率，X为在n次贝努里试验中出现成功的次数，C(X,n)表示在n次试验中出现X的各种组合情况，在此称为二项系数（binomialcoefficient）。内容简介二项分布，伯努里分布：进行一系列试验,如果1. 在每次试验中只有两种可能的结果,而且是互相对立的；2. 每次实验是独立的，与其它各次试验结果无关；3. 结果事件发生的概率在整个系列试验中保持不变，则这一系列试验称为伯努力试验。在这试验中，事件发生的次数为一随机事件,它服从二次分布。二项分布可以用于可靠性试验。可靠性试验常常是投入n个相同的式样进行试验T小时，而只允许k个式样失败，应用二项分布可以得到
26、通过试验的概率。一个事件必然出现，就说它100%要出现。100%=1，所以100%出现的含义就是出现的概率P=1。即必然事件的出现概率为1。若掷一枚硬币，正面向上的结果的概率为0.5。反面向上的结果的概率也是0.5。则出现正面向上事件或者反面向上事件的概率就是0.5+0.5=1，即二者必居其一。若掷两次硬币，由独立事件的概率乘法定理那么两次都是正面（反面）向上的概率是0.50.5=0.25。另外第一个是正第二个是反的出现概率也是0.50.5=0.25。同理第一个反第二个正的出现概率也是0.50.5=0.25。于是一正一反的概率是前面两个情况的和，即0.25+0.25=20.25=0.5。它们的
27、合计值仍然是1。binopdf二项分布的概率密度函数语法格式Y = binopdf(X,N,P)函数功能Y = binopdf(X,N,P) 计算X中每个X(i)的概率密度函数，其中，N中对应的N(i)为试验数，P中对应的P(i)为每次试验成功的概率。Y, N, 和 P 的大小类型相同，可以是向量、矩阵或多维数组。输入的标量将扩展成一个数组，使其大小类型与其它输入相一致。N中的值必须是正整数，0 P 1 。对于给出的x和两个参数n和p，二项分布概率密度函数为其中 q = 1 p。 y为n次独立试验中成功x次的概率，其中，每次试验成功的概率为p。指示器函数 I(0,1,.,n)(x) 确保x 取
28、值为 0, 1, ., n。 示例一个质量检查技术员一天能测试200块电路板。假设有 2% 的电路板有缺陷，该技术员在一天的测试中没有发现有缺陷的电路板的概率是多少？binopdf(0,200,0.02)ans =0.0176质量检查技术员一天中最大可能检测出有缺陷的电路板是多少块？defects=0:200;y = binopdf(defects,200,.02);x,i=max(y);defects(i) ans =4相关函数binocdf | binofit | binoinv | binornd | binostat | pdfdlmread将以ASCII码分隔的数值数据文件读入到矩阵
29、中语法格式M = dlmread(filename)M = dlmread(filename, delimiter)M = dlmread(filename, delimiter, R, C)M = dlmread(filename, delimiter, range)函数功能M = dlmread(filename)读取有分隔符的ASCII数值数据文件filename，并把数据给到矩阵M中。文件名(filename)以字符串形式用单引号括起来。 dlmread 从文件的格式中推断分隔符。M = dlmread(filename, delimiter)指定分隔符delimiter。使用 t 代
30、表制表符 tab 作为分隔。M = dlmread(filename, delimiter, R, C)R和C指定了数据读取数据，其左上角的值位于文件中的第R行，第C列。R和C从0开始，所以R=0, C=0 指向文件中的第一个值。M = dlmread(filename, delimiter, range)读取由range = R1 C1 R2 C2指定的区域块， (R1,C1)是左上角，(R2,C2)是右下角。也可以使用电子表格表示法来指定，如range = A1.B7。备注 All data in the input file must be numeric. dlmread does n
31、ot read files that contain nonnumeric data, even if the specified rows and columns contain only numeric data. 若没有指定分隔符，当从文件格式中推断分隔符时，连续的空格符当作一个分隔符对待。若指定了分隔符，则重复的分隔符将分别作为单独的分隔符对待。 If you want to specify an R, C, or range input, but not a delimiter, set the delimiter argument to the empty string, (two
32、 consecutive single quotes with no spaces in between, ). For example,M = dlmread(myfile.dat, , 5, 2)In this case, dlmread treats repeated white spaces as a single delimiter. Dlmread将用0填充没有边界的区域。有多行的数据文件，若以非空格分隔符结束，例如分号，则导入后会多产生全0的一列在最后。 Dlmread在导入任何复数时，将其作为一个整体导入到一个复数单元中。下面是有效的复数格式：i|jExample: 5.7-3
33、.1ii|jExample: -7j嵌入了空格的复数是不正确的格式，空格将被认为是分隔符。示例例1Export a 5-by-8 test matrix M to a file, and read it with dlmread, first with no arguments other than the filename:M = gallery(integerdata, 100, 5 8, 0); dlmwrite(myfile.txt, M, delimiter, t)dlmread(myfile.txt)ans = 96 77 62 41 6 21 2 42 24 46 80 94 3
34、6 20 75 85 61 2 93 92 82 61 45 53 49 83 74 42 1 28 94 21 90 45 18 90 14 20 47 68Now read a portion of the matrix by specifying the row and column of the upper left corner:dlmread(myfile.txt, t, 2, 3)ans = 92 82 61 45 53 42 1 28 94 21 90 14 20 47 68This time, read a different part of the matrix using
35、 a range specifier:dlmread(myfile.txt, t, C1.G4)ans = 62 41 6 21 2 80 94 36 20 75 93 92 82 61 45 74 42 1 28 94例2Export matrix M to a file, and then append an additional matrix to the file that is offset one row below the first:M = magic(3);dlmwrite(myfile.txt, M*5 M/5, )dlmwrite(myfile.txt, M/3, -ap
36、pend, . roffset, 1, delimiter, )type myfile.txt40 5 30 1.6 0.2 1.215 25 35 0.6 1 1.420 45 10 0.8 1.8 0.4 2.6667 0.33333 21 1.6667 2.33331.3333 3 0.66667When dlmread imports these two matrices from the file, it pads the smaller matrix with zeros:dlmread(myfile.txt) 40.0000 5.0000 30.0000 1.6000 0.200
37、0 1.2000 15.0000 25.0000 35.0000 0.6000 1.0000 1.4000 20.0000 45.0000 10.0000 0.8000 1.8000 0.4000 2.6667 0.3333 2.0000 0 0 0 1.0000 1.6667 2.3333 0 0 0 1.3333 3.0000 0.6667 0 0 0替代作为dlmread的替代，可使用导入向导。选择菜单 “File | Import Data” 激活导入向导。相关函数dlmwrite | textscanfprintf写数据到文本文件或输出到命令窗口语法格式fprintf(fileID,
38、 format, A, .)fprintf(format, A, .)count = fprintf(.)函数功能fprintf(fileID, format, A, .) fileID为文件句柄，指定要写入的文件；format是用来控制所写数据格式的格式符；A是用来存放数据的矩阵。applies the format to all elements of array A and any additional array arguments in column order, and writes the data to a text file. fprintf uses the encodin
39、g scheme specified in the call to fopen.fprintf(format, A, .) 将A中的数据按格式format输出到命令窗口。count = fprintf(.) 返回fprintf写的字节数。输入参数fileID取值为下面中的一种： 从文件打开时所得到的文件句柄，为一个整数值。 1 为标准输出设备（屏幕）。 2 for standard error.fileID缺省时，输出到屏幕，即1为默认值。format单引号括起来的是 “转换控制字符串 ”，指定输出的格式。由下面给出的各部分的组合而成： 百分号后紧跟的格式字符，如%s 。 输出区域的宽度，精度
40、，和其它选项（附加格式说明符）。 要原样输出的文本（普通字符）。 转义字符，包括：单引号% 百分号反斜线a响铃b退格符fForm feedn换行r回车t水平制表符v垂直制表符xN十六进制数, NN八进制数, N转换字符和可选的操作符按以下顺序排列（包括空格）：下表列出可用的转换字符和限定类型字符。数据类型转换说 明Integer, signed%d 或 %i十进制整数%ld 或 %li64-bit base 10 values%hd 或 %hi16-bit base 10 valuesInteger, unsigned%u十进制整数%o八进制整数(octal)%x十六进制整数 (hexadec
41、imal), 使用字母 af%X十六进制整数 (hexadecimal), 使用字母AF%lu%lo%lx or %lX64-bit values, 十、八、十六进制整数%hu%ho%hx or %hX16-bit values, 十、八、十六进制整数Floating-point number%f定点表示法%e指数表示法，如 3.141593e+00%E同 %e，但用E，如 3.141593E+00%g%e 或 %f的紧凑格式，没有尾部的0%G%E 或 %f的紧凑格式，没有尾部的0%bx or %bX%bo%bu双精度的十六、八、十进制数值如：%bx prints pi as 400921fb
42、54442d18%tx or %tX%to%tu单精度的十六、八、十进制数值如：%tx prints pi as 40490fdbCharacters%c单个字符%s字符串附加格式说明符，包括： 字段宽度Minimum number of characters to print. Can be a number, or an asterisk (*) to refer to an argument in the input list. For example, the input list (%12d, intmax) is equivalent to (%*d, 12, intmax). 精
43、度For %f, %e, or %E: Number of digits to the right of the decimal point.Example: %6.4f prints pi as 3.1416For %g or %G Number of significant digits.Example: %6.4g prints pi as 3.142Can be a number, or an asterisk (*) to refer to an argument in the input list. For example, the input list (%6.4f, pi) i
44、s equivalent to (%*.*f, 6, 4, pi). FlagsActionFlagExampleLeft-justify.%-5.2fPrint sign character (+ or ).+%+5.2fInsert a space before the value. % 5.2fPad with zeros.0%05.2fModify selected numeric conversions: For %o, %x, or %X, print 0, 0x, or 0X prefix. For %f, %e, or %E, print decimal point even wh
