【導(dǎo)讀】能夠證明某些模型或理論是正確。時(shí)間序列的例子。的預(yù)測(cè)，如何點(diǎn)圖,……通常是等間距的:. 確定的或隨機(jī)的(精確預(yù)測(cè)是不可。人們對(duì)統(tǒng)計(jì)數(shù)據(jù)往往可以根據(jù)其特點(diǎn)從兩。個(gè)方面來(lái)切入，以簡(jiǎn)化分析過(guò)程。對(duì)象的觀測(cè)值組成的數(shù)據(jù)。另一個(gè)稱為時(shí)間序列，也就是。前面討論的模型多是和橫截面數(shù)據(jù)有關(guān)。論更加復(fù)雜的包含這兩方面的數(shù)據(jù)。

　　

【正文】 w(y) A=matrix(1,2,1) A is a 2x1 matrix of 1s same for all t mu0= Sigma0= .01 Phi=1 initpar=c(.1,.1,.1,0) initial parameter values in order, para[1] to para[4] Estimation view help(optim) for details here... the iteration number is written to the screen but you have to manually scroll down to see it est=optim(initpar,Linn,NULL,method=BFGS,hessian=TRUE,control=list(trace=1,REPORT=1)) stderr=sqrt(diag(solve(est$hessian))) display summary of estimation estimate=est$par u=cbind(estimate,stderr) rownames(u)=c(sigw,cR11, cR22, cR12) u Smooth first set parameters to their final estimates cQ=est$par[1] cR1=est$par[2] 11 element of chol(R) cR2=est$par[3] 22 element of chol(R) cR12=est$par[4] 12 element of chol(R) cR=matrix(c(cR1,0,cR12,cR2),2) ks=Ksmooth0(num,y,A,mu0,Sigma0,Phi,cQ,cR) Plot ((ks$xs),lwd=2,ylim=c(,.4), ylab=Temp Deviations) lines(y1,col=blue,lty=dashed) color helps here lines(y2,col=red, lty=dashed) estimate stderr sigw cR11 cR22 cR12 !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! assumes ss0 has been sourced !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! generate data (999) num=100 N=num+1 need 101 x39。s generate x(t)=.8x(t1)+w(t), t=1,...,100 and x(0) from the stationary distribution: x = (n=N, list(ar = .8, sd=1)) here you have x0 to x100 v = rnorm(num,0,1) obs noise y=ts(x[1]+v) observations y[1] ,..., y[100] initial estimates (method discussed in the text) u=(y,lag(y,1),lag(y,2)) varu=var(u)。 coru=cor(u) phi=coru[1,3]/coru[1,2] q = (1phi^2)*varu[1,2]/phi r = varu[1,1] q/(1phi^2)。 initpar=c(phi,sqrt(q),sqrt(r)) initpar view the initial estimates function to evaluate innovations likelihood Linn=function(para){ phi=para[1] sigw=para[2] this is the standard dev sigv=para[3] this is the standard dev mu0=0 Sigma0=(sigw^2)/(1phi^2) Sigma0[Sigma00] =0 make sure Sigma0 is never negative kf = Kfilter0(num,y,1,mu0,Sigma0,phi,sigw,sigv) run filter under present parameters return(kf$like) return log likelihood } do the estimation view help(optim) for details here... the iteration number is written to the screen but you have to manually scroll down to see it est=optim(initpar, Linn, NULL, method = BFGS, hessian = TRUE, control=list(trace=1,REPORT=1)) stderr=sqrt(diag(solve(est$hessian))) standard errors display summary of estimation estimate=est$par u=cbind(estimate,stderr) rownames(u)=c(phi,sigw,sigv) u Run estimate stderr phi sigw sigv Global Warming T i m eTemp Deviations0 20 40 60 80 1000.60.40.20.00.20.4 via BFGS [167。 ]: !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! assumes ss0 has been sourced !!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!!! Read Data amp。 Make Measurement Matrix y=ts(scan(c:/xzwu/2020/berkeley/ts/shumway/mydata/),start=1960,freq=4) num=length(y) A=cbind(1,1,0,0) Function to Calculate Likelihood Linn=function(para){ phi=para[1] Phi=diag(0,4)。 Phi is 4x4 but only one element is a parameter Phi[1,1]=phi。 Phi[2,]=c(0,1,1,1) Phi[3,] = c(0,1,0,0)。 Phi[4,] = c(0,0,1,0) cQ1=para[2] sqrt q11 cQ2=para[3] sqrt q22 cQ=diag(0,4)。 cQ[1,1]=cQ1。 cQ[2,2]=cQ2 cR=para[4] sqrt r11 kf=Kfilter0(num,y,A,mu0,Sigma0,Phi,cQ,cR) return(kf$like) } Initial Parameters mu0=c(.7,0,0,0) Sigma0= diag(.04,4) initpar=c(,.1,.1,.5) initial parameters for Phi[1,1], the 2 Q39。s and R Estimation the iteration number is printed to the screen but you have to manually scroll to see it est=optim(initpar,Linn,NULL,method=BFGS,hessian=TRUE,control=list(trace=1,REPORT=1)) stderr=sqrt(diag(solve(est$hessian))) display summary of estimation estimate=est$par u=cbind(estimate,stderr) rownames(u)=c(Phi11,sigw1,sigw2,sigv) u Smooth phi=est$par[1] Phi=diag(0,4)。 Phi[1,1]=phi。 Phi[2,]=c(0,1,1,1) Phi[3,] = c(0,1,0,0)。 Phi[4,] = c(0,0,1,0) cQ1=est$par[2] cQ2=est$par[3] cQ=diag(1,4)。 cQ[1,1]=cQ1。 cQ[2,2]=cQ2 note lower diag is 2x2 ident for inversions (as a device, they39。re not used) cR=est$par[4] ks=Ksmooth0(num,y,A,mu0,Sigma0,Phi,cQ,cR) Plot Tsm=ts(ks$xs[1,],start=1960,freq=4) Ssm=ts(ks$xs[2,],start=1960,freq=4) p1=2*sqrt(ks$Ps[1,1,]) p2=2*sqrt(ks$Ps[2,2,]) par(mfrow=c(3,1)) plot(Tsm, main=Trend Component, ylab=Trend) lines(Tsm+p1,lty=dashed, col=blue) lines(Tsmp1,lty=dashed, col=blue) plot(Ssm, main=Seasonal Component, ylim=c(5,4), ylab=Season ) lines(Ssm+p2,lty=dashed, col=blue) lines(Ssmp2,lty=dashed, col=blue) plot(y, type=p, main=Data (points) and Trend+Season (line)) lines(Tsm+Ssm) estimate stderr Phi11 sigw1 sigw2 sigv Tre nd C om po ne ntT i m eTrend1960 1965 1970 1975 198051015S e a s on a l C om po ne ntT i m eSeason1960 1965 1970 1975 19804024D a t a ( po i nt s ) a nd Tre nd + S e a s on ( l i ne )y1960 1965 1970 1975 1980051015dse formulas or（ 更一般 where n_t is the system noise, Q, the system noise matrix, and R the output (measurement) noise matrix.) ARMA [including ARIMA, VAR: B(L)=I] State Space(A linear timeinvariant state space representation in innovations form i) output input Lag operator 不可觀測(cè)的狀態(tài)向量 注意：看 以及