60 likes | 131 Views
Lecture: EPS 236, 17 Sep 2012. This lecture: Autoregressive time series; variance; relationship to a Markov process. [Examples: 1 st order cases] [autoregressive_var.r] Coming next: Fitting lines to data—ordinary least squares
E N D
Lecture: EPS 236, 17 Sep 2012 This lecture: • Autoregressive time series; variance; relationship to a Markov process. [Examples: 1st order cases][autoregressive_var.r] Coming next: • Fitting lines to data—ordinary least squares • Why ordinary least squares is often incorrect as applied to our data sets [ols.r] 1
Discrete 1st order Markov process Note the larger variance of the autocorrelated time series, also the prevalence of the decorrelation time (10 yr). θ=.904 N=1000 tt=1:N k=.10 theta = exp(-k) a=rnorm(N) quartz() # X11( ) plot(tt,a) x=0 for(i in 1:N) {x=c(x,x[i]*theta + a[i])} #var(x)=5.5 ; var(a)=.92; 1/(1- θ2)=5.5 X11( ) plot(c(0,tt),x,type="l") X11( ) plot(c(0,tt),x,type="l",xlim=c(100,200), ylim=c(-5,5),err=-1) axis(side=3,col="red",at=seq(100,200,10), tck=.1,labels=F) abline(v=seq(100,200,10),lty=2,col="red")
X Prob x > X p(x)=2/(σ√π)exp(-x2/σ2) prob X< x' 2/(σ√π)∫exp(-t2/σ2)dt Prob x < X units of σ
tr1 = function(x,p1){ if(x > p1) ans =1 else ans = 0 ans} #symmetrical 2-state markov chain; eqm, => equal pops. # transition probability matrix -- markov2 = function(nnn,p1=.3){ #2 state markov process dum = runif(nnn+1) ##not rnorm state = 1 + as.numeric(rnorm(1)>0) # random initial state state.cum = NULL for(ii in dum) { if(state == 1) { state = state + tr1(ii,p1) } if(state == 2) { state = state - 1 + tr1(ii,p1) } state.cum = c(state.cum, state) } state.cum1= state.cum[(1:nnn)+ 1] dx = state.cum1-state.cum[1:nnn] # detect if we have a transition pp = c(sum(state.cum[1:nnn]==1&dx==0),sum(dx==1),sum(dx==-1), sum(state.cum[1:nnn]==2&dx==0)) pz = c(pp[1]/(pp[1]+pp[2]),pp[2]/(pp[1]+pp[2]),pp[3]/(pp[3]+ pp[4]),pp[4]/(pp[3]+pp[4])) ans = c(p1,nnn,pp,pz) names(ans) = c("prob","n","n11","n12","n21","n22","p11","p12","p21","p22") ans } # markov.res = NULL for(nn in c(rep(10,3000), rep(25,1000), rep(50,500), rep(250,100), rep(1000,60), rep(5000,40))) { markov.res = rbind(markov.res, markov2(nn)) } zum = tapply(markov.res[,"p11"]-.3, markov.res[,"n"],var,na.rm=T) p11 1 p12 p21 2 p22