######################################################################################################### swaveqexMerge <- function(cdatin,qwstnum,ddstnum,yrbeg,yrend,getdd="WD") { # # merges pesticide data with daily discharge data and prepares # input file for swaveqexFit # # cdatin is a data frame with the pesticide concentration data # - the first column should be the station number (character) # - the second column should be the date, in "yyyy-mm-dd" format (character) # - the concentration value should be a column named "final_value" (numeric) # - the remark should be in a column names "final_remark" (character) # # qwstnum is the station number to analyze # ddstnum is the station number for daily discharge (usually the same as qwstnum) # yrbeg and yrend are the beginning and ending calendar years for analysis (numeric) # # getdd specifies where to get the daily discharge data # getdd="WD" (the default) uses the waterData package to get streamflow # data from the USGS web server # getdd="File" reads daily discharge from the tab-delimited text file # called "dd_ddstnum.txt" in the default directory # # # Requires library waterData # # select pesticide data for specified site and # download daily discharge data using importDVs from the waterData library cdatxx <- cdatin[cdatin[,1]==qwstnum,] cdatxx <- data.frame(cdatxx[,2],cdatxx[,"final_remark"],cdatxx[,"final_value"]) names(cdatxx) <- c("date","crem","cval") yrxx <- as.numeric(substring(as.character(cdatxx[,1]),1,4)) if(min(yrxx)>yrbeg) {yrbeg <- min(yrxx) } if(max(yrxx)=yrbeg & yr<=yrend,] qdate <- as.character(qin[,1]) qval <- qin[,2] qval[qval<1] <- 1 yr <- as.numeric(substring(qdate,1,4)) mo <- as.numeric(substring(qdate,6,7)) da <- as.numeric(substring(qdate,9,10)) qday <- rep(NA,length(qdate)) cumd <- c(0,31,59,90,120,151,181,212,243,273,304,334) mday <- c(31,28,31,30,31,30,31,31,30,31,30,31) chmo <- c("01","02","03","04","05","06","07","08","09","10","11","12") chday <- c(paste("0",as.character(1:9),sep=""),as.character(10:31)) for (j in 1:12) { qday[mo==j] <- (yr[mo==j]-yrbeg)*365 + cumd[j]+da[mo==j] } qdaily <- rep(NA,365*(yrend-yrbeg+1)) qdaily[qday] <- qval ntot <- 365*(yrend-yrbeg+1) qdate <- rep(" ",ntot) for (i in yrbeg:yrend) { for (j in 1:12) { for (k in 1:mday[j]) { qdate[(i-yrbeg)*365+cumd[j]+k] <- paste(as.character(i),"-",chmo[j],"-",chday[k],sep="") } } } # compute flow anomalies qmean <- mean(log10(qdaily),na.rm=T) qmc <- log10(qdaily)- qmean qanom30 <- rep(0,nall) for (i in 30:nall) {qanom30[i] <- mean(qmc[(i-29):i])} qanom30[1:29] <- qanom30[30] qanom1 <- qmc-qanom30 datxxx <- data.frame(qdate,1:nall,rep(1:365,rlen),qanom30,qanom1) datxxx[,1] <- as.character(datxxx[,1]) names(datxxx) <- c('date','jday','sday','qanom30','qanom1') # # prepare concentration data # chdate2 <- as.character(cdatxx[,1]) yr2 <- as.numeric(substring(chdate2,1,4)) mo2 <- as.numeric(substring(chdate2,6,7)) da2 <- as.numeric(substring(chdate2,9,10)) leapdays <- (mo2==2 & da2==29) cdatxx <- cdatxx[!leapdays & yr2>=yrbeg & yr2<=yrend,] chdate2 <- as.character(cdatxx[,1]) yr2 <- as.numeric(substring(chdate2,1,4)) mo2 <- as.numeric(substring(chdate2,6,7)) da2 <- as.numeric(substring(chdate2,9,10)) jday2 <- (yr2-yrbeg)*365+ cumd[mo2] + da2 sday2 <- cumd[mo2] + da2 #### thin to at least three days between samples nx2 <- length(jday2) pkthin <- rep(F,nx2) pkthin[1] <- T jdaycum <- jday2[1] for (ii in 2:nx2) { if ((jday2[ii]-jdaycum)>2) { pkthin[ii] <- T jdaycum <- jday2[ii] } } cdatxx <- cdatxx[pkthin,] jday2 <- jday2[pkthin] sday2 <- sday2[pkthin] # remove values with missing streamflow anomalies qchkna <- is.na(qanom30[jday2]) | is.na(qanom1[jday2]) cdatxx <- cdatxx[!qchkna,] jday2 <- jday2[!qchkna] sday2 <- sday2[!qchkna] # # create dataframes for input to swaveqexFit # nobs <- length(jday2) icen2xx <- as.character(cdatxx[,2]) icen2 <- rep("u",nobs) icen2[icen2xx=="<"] <- "c" cdatxx[,2] <- as.character(icen2) datyyy <- data.frame(cdatxx[,1],jday2,sday2,log10(cdatxx[,3]),cdatxx[,2],qanom30[jday2],qanom1[jday2]) datyyy[,1] <- as.character(datyyy[,1]) datyyy[,5] <- as.character(datyyy[,5]) names(datyyy) <- c('date','jday','sday',"lgconc","icen",'qanom30','qanom1') cenyr <- nall/365/2 icenall <- rep("u",nall) icenall[jday2[datyyy[,"icen"]=="c"]] <- "c" datsub <- data.frame(datyyy[,1:5],datyyy[,4],rep(1,nobs),rep(NA,nobs),datyyy[,6:7],datyyy[,2]/365-cenyr) datall <- data.frame(datxxx[,1:3],rep(NA,nall),icenall,rep(NA,nall),rep(1,nall),rep(NA,nall),datxxx[,4:5],datxxx[,2]/365-cenyr) datsub[,1] <- as.character(datsub[,1]) datsub[,5] <- as.character(datsub[,5]) datall[,1] <- as.character(datall[,1]) datall[,5] <- as.character(datall[,5]) names(datsub) <- c('date','jday','sday','yobs','icen','yimp','int','swave','qanom30','qanom1','tnd') names(datall) <- c('date','jday','sday','yobs','icen','yimp','int','swave','qanom30','qanom1','tnd') # produce rough plot of observations and flow anomalies xxx <- yrbeg+datsub[,2]/365 yyy <- datsub[,4] icen <- as.character(datsub[,5]) plot(xxx,yyy,xlim=c(yrbeg,yrend+1),type="p",pch=" ") mtext(side=3,line=1,cex=.8,"Observations") for (i in yrbeg:yrend) {abline(v=i) } points(xxx[icen=="u"],yyy[icen=="u"],pch=16,cex=.7) points(xxx[icen=="c"],yyy[icen=="c"],pch=1,cex=.7,col=2) nucen <- sum(icen=="u") ntot <- length(icen) pucen <- 100*round(nucen/ntot,2) mtext(side=3,line=0,paste("numobs=",as.character(ntot)," nuncen=",as.character(nucen),"prucen=",as.character(pucen))) xxx <- yrbeg+datall[,2]/365 yyy <- datall[,9] plot(xxx,yyy,xlim=c(yrbeg,yrend+1),type="l") mtext(side=3,line=1,cex=.8,"Mid-term Flow Anomaly") for (i in yrbeg:yrend) {abline(v=i) } yyy <- datall[,10] plot(xxx,yyy,xlim=c(yrbeg,yrend+1),type="l") mtext(side=3,line=1,cex=.8,"Short-term Flow Anomaly") for (i in yrbeg:yrend) {abline(v=i) } sdqamt <- sd(datall[,"qanom30"],na.rm=T) sdqast <- sd(datall[,"qanom1"],na.rm=T) saveout <- list(qwstnum,c(yrbeg,yrend),datsub,datall,qdaily,c(ntot,nucen,pucen),c(sdqamt,sdqast)) } ######################################################################################################### swaveqexFit <- function(datin,outfolder,ncs=50) { # # Fits the seawaveqex model using input data prepared by swaveqexMerge # Produces plot file with diagnostic plots and text file with conditional simulations # # datin is a list produced by swaveqexMerge # outfolder is a character variable specifying the folder to put the output files # For example, "outAtrazine\\" puts the files in a folder named outAtrazine in # the default directory. The folder needs to be created beforehand. # qwstnum <- datin[[1]] yrbeg <- datin[[2]][1] yrend <- datin[[2]][2] datsub <- datin[[3]] datall <- datin[[4]] qdaily <- datin[[5]] qdate <- as.character(datall[,1]) pdf(file=paste(outfolder,"Plots",qwstnum,".pdf",sep=""),height=6,width=10) par(omi=c(.2,.2,.5,.2)) rlen <- yrend-yrbeg+1 nall <- rlen*365 svcounts <- datin[[6]] svsdflowanoms <- datin[[7]] # # Seasonal wave selection and parameter estimation # svstuff <- swaveqexPESTpdo(datsub) # Generate conditional simulations # ngen is the number of conditional simulations # parxxx <- svstuff[1:15] spmodcls <- svstuff[16] svpvals <- svstuff[17:20] svstuff <- svstuff[1:15] ngen <- ncs if(ngen>250) {ngen <- 250} tmpout <- swaveqexCSIM(ngen,parxxx,datsub,datall,spmodcls) estannmax <- tmpout[[1]] datsub <- data.frame(datsub,tmpout[[2]]) datall <- data.frame(datall,tmpout[[3]]) svcondsim <- round(10^tmpout[[4]],5) for (jjj in 1:ngen) { pickxxx <- svcondsim[,jjj] < 0.00001 svcondsim[pickxxx,jjj] <- 0.00001 pickxxx <- is.na(svcondsim[,jjj]) svcondsim[pickxxx,jjj] <- (-0.009) } svcenout <- rep(" ",nall) svcenout[datsub[,"jday"]] <- datsub[,"icen"] svcenout[svcenout=="c"] <- "<" svcenout[svcenout=="u"] <- "o" concorig <- (round(10^datsub[,"yobs"],5)) svconcout <- rep(" ",nall) svconcout[datsub[,"jday"]] <- as.character(concorig) qdate <- as.character(datall[,1]) qyear <- as.numeric(substring(qdate,1,4)) qsday <- datall[,3] qout <- round(qdaily,1) qout[is.na(qout)] <- (-9.0) muout <- round(10^datall[,"mu"],5) mucondout <- round(10^datall[,"condmean"],5) muout[muout<0.00001] <- 0.00001 mucondout[mucondout<0.00001] <- 0.00001 muout[is.na(muout)] <- (-0.009) mucondout[is.na(mucondout)] <- (-0.009) svcondsim <- 1000*svcondsim muout <- 1000*muout mucondout <- 1000*mucondout # # save conditional simulations in a file # savecsim <- data.frame(qdate,qyear,qsday,qout,svconcout,svcenout,muout,mucondout,svcondsim[,1:min(ngen,100)],stringsAsFactors=F) for (jjj in 7:(9+min(ngen,100)-1) ) { savecsim[,jjj] <- format(savecsim[,jjj],scientific=F) } # savecsim[,1] <- as.character(savecsim[,1]) # savecsim[,5] <- as.character(savecsim[,5]) # savecsim[,6] <- as.character(savecsim[,6]) names(savecsim) <- c("date","year","jday","qobs","cobs","crem","estreg","estcmu",paste("csim",as.character(1:min(ngen,100)),sep="")) write.table(format(savecsim),file=paste(outfolder,"CSIMS",qwstnum,".txt",sep=""),row.names=F,sep="\t",quote=F) # # produce plots # hlife <- parxxx[10] ; cmax <- parxxx[8] ; modn <- parxxx[9] swave <- compwaveconvxx(cmax,modn,hlife,spmodcls) sdum <- (1:365)/366 ipk <- floor(360*sdum) ipk[ipk==0] <- 1 swave <- swave[ipk] swaveobs <- swave[datsub[,"sday"]] swaveall <- swave[datall[,"sday"]] datsub[,"swave"] <- swaveobs datall[,"swave"] <- swaveall annmax <- log10(estannmax) condsim <- datall[,"condsim"] condmean <- datall[,"condmean"] yobs <- datsub[,"yobs"] condsimobs <- datsub[,"condsim"] yymax <- ceiling(max(c(annmax,condsim,condmean)+0.1,na.rm=T)) yymin <- floor(min(c(annmax,condsim,condmean),na.rm=T)) if(yymin< (-6)) {yymin <- (-6) } jdayobs <- datsub[,"jday"] ; dyrobs <- yrbeg+(jdayobs)/365 jdayall <- datall[,"jday"] ; dyrall <- yrbeg+(jdayall)/365 sdayobs <- datsub[,"sday"] ; sdayall <- datsub[,"sday"] muobs <- datsub[,"mu"] muall <- datall[,"mu"] omgobs <- datsub[,"omg"] omgall <- datall[,"omg"] yhatobs <- datsub[,"yhat"]; yhatall <- datall[,"yhat"] icen <- datsub[,"icen"]=="c" par(mai=c(.8,1,.2,.2)) plot(dyrall,condsim,type="p",pch=" ",xlim=c(yrbeg-0.2,yrend+1.2),ylim=c(yymin,yymax),xaxt="n",yaxt="n",ylab=" ",xlab=" ",yaxs="i",xaxs="i") lines(dyrall,condsim,type="l",col="dark grey") for (j in 0:(rlen)) {abline(v=yrbeg+j,col=5,lty=2)} for (j in 0:(rlen-1)) { lines(c(yrbeg+j,yrbeg+j+1),c(annmax[j+1],annmax[j+1]),col=5,lwd=2) mtext(side=1,line=0.5,cex=0.7,at=yrbeg+j+0.5,as.character(yrbeg+j)) } points(dyrobs[!icen],condsim[jdayobs[!icen]],pch=16,col="red",cex=0.6) points(dyrobs[icen],condsim[jdayobs[icen]],pch=1,col="red",cex=0.7) if (sum(icen)>0) { abline(h=min(yobs[icen]),col="red") } for (j in yymin:yymax) { mtext(side=2,line=0.5,adj=1,cex=0.7,at=j,as.character(10^j),las=1) lines(c(yrbeg-0.2,yrbeg-0.2+0.015*rlen),c(j,j)) ; lines(c(yrend+1.2-0.015*rlen,yrend+1.2),c(j,j)) } yydel <- (yymax-yymin)*0.025 if (sum(icen)==0) {polygon(c(yrbeg+.31*rlen,yrbeg+.59*rlen,yrbeg+.59*rlen,yrbeg+.31*rlen), c(yymax-6*yydel,yymax-6*yydel,yymax-yydel,yymax-yydel),col="white")} if (sum(icen)>0) {polygon(c(yrbeg+.31*rlen,yrbeg+.59*rlen,yrbeg+.59*rlen,yrbeg+.31*rlen), c(yymax-8*yydel,yymax-8*yydel,yymax-yydel,yymax-yydel),col="white")} text(yrbeg+.45*rlen,yymax-2*yydel,cex=.6,"EXPLANATION") lines(c(yrbeg+.32*rlen,yrbeg+.38*rlen),rep(yymax-3*yydel,2),col="dark grey",lwd=1.5) text(yrbeg+.4*rlen,yymax-3*yydel,cex=.6,adj=0,"Conditional trace") lines(c(yrbeg+.32*rlen,yrbeg+.38*rlen),rep(yymax-4*yydel,2),col=5,lwd=2) text(yrbeg+.4*rlen,yymax-4*yydel,cex=.6,adj=0,"Estimated annual maximum") points(yrbeg+0.35*rlen,yymax-5*yydel,pch=16,col="red",cex=0.6) text(yrbeg+.4*rlen,yymax-5*yydel,cex=.6,adj=0,"Observed concentration") if (sum(icen)>0) { points(yrbeg+0.35*rlen,yymax-6*yydel,pch=1,col="red",cex=0.7) text(yrbeg+.4*rlen,yymax-6*yydel,cex=.6,adj=0,"Simulated censored concentration") lines(c(yrbeg+.32*rlen,yrbeg+.38*rlen),rep(yymax-7*yydel,2),col="red") text(yrbeg+.4*rlen,yymax-7*yydel,cex=.6,adj=0,"Censoring limit") } mtext(side=1,line=2,cex=0.8,"Calendar year") mtext(side=2,line=3,cex=0.8,"Concentration, in micrograms per liter") # mtext(side=3,outer=T,line=0,cex=1.2,paste(outfolder,qwstnum)) par(mai=c(.8,1.5,.5,.7)) yy <- condsimobs-yhatobs plot(sdayobs,yy,type="p",pch=" ",xlim=c(-15,380),ylim=c(yymin,yymax),xaxt="n",yaxt="n",ylab=" ",xlab=" ",yaxs="i",xaxs="i") points(sdayobs[!icen],yy[!icen],pch=16,col="red",cex=0.6) points(sdayobs[icen],yy[icen],pch=1,col="red",cex=0.7) lines(1:365,svstuff[1]+svstuff[2]*swaveall[1:365],lwd=2) lines(1:365,svstuff[1]+svstuff[2]*swaveall[1:365]+2*omgall[1:365]*svstuff[14],lwd=2,lty=4) lines(1:365,svstuff[1]+svstuff[2]*swaveall[1:365]-2*omgall[1:365]*svstuff[14],lwd=2,lty=4) for (j in yymin:yymax) { mtext(side=2,line=0.5,adj=1,cex=0.7,at=j,as.character(10^j),las=1) lines(c(-15,-10),c(j,j)) ; lines(c(375,380),c(j,j)) } mocum <- c(0,31,59,90,120,151,181,212,243,273,304,334) for (j in c(mocum,365)) { abline(v=j,col=5,lty=2) } mtext(side=1,line=0.5,cex=0.7,at=mocum+15,c("Jan.","Feb.","Mar.","Apr.","May","June","July","Aug.","Sept.","Oct.","Nov.","Dec.")) if (sum(icen)==0) {polygon(c(.61*365,370,370,.61*365), c(yymax-7*yydel,yymax-7*yydel,yymax-yydel,yymax-yydel),col="white")} if (sum(icen)>0) {polygon(c(.61*365,370,370,.61*365), c(yymax-8*yydel,yymax-8*yydel,yymax-yydel,yymax-yydel),col="white")} text(.75*365,yymax-2*yydel,cex=.6,"EXPLANATION") lines(c(.62*365,.68*365),rep(yymax-3*yydel,2),lty=1,lwd=2) text(.7*365,yymax-3*yydel,cex=.6,adj=0,"Seasonal wave") lines(c(.62*365,.68*365),rep(yymax-4*yydel,2),lty=4,lwd=2) text(.7*365,yymax-4*yydel,cex=.6,adj=0,"Seasonal wave plus and minus two times SSD") points(0.65*365,yymax-5*yydel,pch=16,col="red",cex=0.6) text(.7*365,yymax-5*yydel,cex=.6,adj=0,"Adjusted concentration") if (sum(icen)==0) { text(.7*365,yymax-6*yydel,cex=.6,adj=0,"SSD, seasonal standard deviation") } if (sum(icen)>0) { points(0.65*365,yymax-6*yydel,pch=1,col="red",cex=0.7) text(.7*365,yymax-6*yydel,cex=.6,adj=0,"Simulated censored adjusted concentration") text(.7*365,yymax-7*yydel,cex=.6,adj=0,"SSD, seasonal standard deviation") } mtext(side=2,line=3,cex=0.8,"Adjusted concentration, in micrograms per liter") # mtext(side=3,outer=T,line=0,cex=1.2,paste(outfolder,qwstnum)) yy <- (condsimobs-muobs)/omgobs/svstuff[14] yy[yy>3] <- 3 yy[yy< (-3)] <- (-3) plot(sdayobs,yy,type="p",pch=" ",xlim=c(-15,380),ylim=c(-3.2,3.2),xaxt="n",yaxt="n",ylab=" ",xlab=" ",yaxs="i",xaxs="i") for (j in (-3):3) { mtext(side=2,line=0.5,adj=1,cex=0.7,at=j,as.character(j),las=1) if (abs(j)<=3) {abline(h=j,col=5,lty=2) } } points(sdayobs[!icen],yy[!icen],pch=16,col="red",cex=0.6) points(sdayobs[icen],yy[icen],pch=1,col="red",cex=0.7) abline(h=0) for (j in c(mocum,365)) { abline(v=j,col=5,lty=2) } mtext(side=1,line=0.5,cex=0.7,at=mocum+15,c("Jan.","Feb.","Mar.","Apr.","May","June","July","Aug.","Sept.","Oct.","Nov.","Dec.")) # if (sum(icen)>0) { # points(0.67*365,2.8,pch=1,col="red",cex=0.7) ; text(.7*365,2.8,cex=.6,adj=0,"Simulated censored residual") # } mtext(side=2,line=3,cex=0.8,"Normalized residual (dimensionless)") # mtext(side=3,outer=T,line=0,cex=1.2,paste(outfolder,qwstnum)) plot(dyrobs,yy,type="p",pch=" ",xlim=c(yrbeg-0.2,yrend+1.2),ylim=c(-3.2,3.2),xaxt="n",yaxt="n",ylab=" ",xlab=" ",yaxs="i",xaxs="i") points(dyrobs[!icen],yy[!icen],pch=16,col="red",cex=0.6) points(dyrobs[icen],yy[icen],pch=1,col="red",cex=0.7) abline(h=0) for (j in (-3):3) { mtext(side=2,line=0.5,adj=1,cex=0.7,at=j,as.character(j),las=1) if (abs(j)<=3) {abline(h=j,col=5,lty=2) } } for (j in 0:(rlen-1)) { abline(v=yrbeg+j+1,col=5,lty=2) mtext(side=1,line=0.5,cex=0.7,at=yrbeg+j+0.5,as.character(yrbeg+j)) abline(v=yrbeg,col=5,lty=2) } # if (sum(icen)>0) { # points(yrbeg+0.67*rlen,2.8,pch=1,col="red",cex=0.7) ; text(yrbeg+.7*rlen,2.8,cex=.6,adj=0,"Simulated censored residual") # } mtext(side=2,line=3,cex=0.8,"Normalized residual (dimensionless)") # mtext(side=3,outer=T,line=0,cex=1.2,paste(outfolder,qwstnum)) ###### correlogram rxx <- yy/sd(yy) txx <- dyrobs nxx <- length(rxx) nxx2 <- nxx*(nxx-1)/2 dtmp <- rep(NA,nxx2) vtmp <- rep(NA,nxx2) cxx <- 0 for(ixx in 1:(nxx-1)) { for(jxx in (ixx+1):nxx) { cxx <- cxx+1 dtmp[cxx] <- txx[jxx]-txx[ixx] vtmp[cxx] <- (rxx[jxx]-rxx[ixx])^2 } } dtmp <- dtmp*365 pcktmp <- dtmp<50 vtmp <- vtmp[pcktmp] dtmp <- dtmp[pcktmp] # plot(dtmp,vtmp) # lines(supsmu(dtmp,vtmp)) evxx <- seq(3.5,45.5,3.5) ; nxx <- length(evxx) ; hwid <- rep(3.5,nxx) evyy <- rep(NA,nxx) for (j in 1:nxx) { evpck <- dtmp>(evxx[j]-hwid[j]) & dtmp<(evxx[j]+hwid[j]) if(sum(evpck)>=5) {evyy[j] <- mean(vtmp[evpck]^0.25)^4/(0.457+0.494/sum(evpck)) } } evyy <- 1 - 0.5*evyy rhoxx <- parxxx[13] ylow <- -0.4 ; yup <- 1.2 ytck <- .02*(yup-ylow) xlow <- 0; xup <- 50; xtck <- .012*(xup-xlow) par(mai=c(0.8,2,.5,1.2)) plot(evxx,evyy,pch=' ',cex=.8,xlim=c(xlow,xup),ylim=c(ylow,yup),xaxs='i',yaxs='i',xaxt='n',yaxt='n',xlab='',ylab='') for (j in seq(xlow,xup,10)) { mtext(side=1,line=.5,at=j,cex=.8,as.character(j)) lines(c(j,j),c(ylow,ylow+ytck),lwd=1); lines(c(j,j),c(yup,yup-ytck),lwd=1) } for (j in seq(ylow,yup,.2)) {mtext(side=2,line=.5,at=j,adj=1,las=1,cex=.8,as.character(j)) lines(y=c(j,j),x=c(xlow,xlow+xtck),lwd=1); lines(y=c(j,j),x=c(xup,xup-xtck),lwd=1) } lines(evxx,evyy,type='p',pch=16,cex=.8,col=1) lines(seq(0,xup,.1),rhoxx^seq(0,xup,.1),lwd=2,col=1) abline(h=0) mtext(side=2,line=3,cex=.8,'Correlation') mtext(side=1,line=2,cex=.8,'Time between observations, in days') text(40,1.1,cex=.7,"EXPLANATION") lines(c(30,33),c(1.05,1.05)) text(34,1.05,adj=0,cex=.7,"Fitted exponential correlation function") points(31.5,1.0,cex=.8,pch=16) text(34,1,adj=0,cex=.7,"Empirical correlogram") polygon(c(29.2,48.6,48.6,29.2),c(.97,.97,1.13,1.13)) dev.off() svstuff[13] <- -1/log(svstuff[13]) svstuff <- round(svstuff,4)[c(1:5,8:10,12:15)] svstuff <- c(svstuff[1:5],spmodcls,svstuff[c(7,8,6,11,9,10,12)]) # estannmax <- round(estannmax,4) # chnames1 <- c("int","cswave","cmtfa","cstfa","ctnd","wmcls","wmodno","hlife","wshft","sigma","alph","cts","neg2LLIK") # chnames2 <- c(paste("cmax",as.character(yrbeg:yrend),sep="")) # svstuff <- data.frame(matrix(c(svstuff,estannmax),nrow=1)) # names(svstuff) <- c(chnames1,chnames2) # row.names(svstuff) <- NULL svstuff <- c(svstuff,svcounts,svpvals,svsdflowanoms) svstuff <- c(svstuff[1:20],yrbeg,yrend,rlen,svstuff[21:22]) svnames <- c("int","cswave","cmtfa","cstfa","ctnd","wmcls","wmodno","hlife","wshft","sigma","alph","cts","n2LLIK", "nobs","nucen","prucen","pswave","pmtfa","pstfa","ptnd","yrbeg","yrend","rlen","sdmtfa","sdstfa") svstuff <- svstuff[c(21,22,23,14,15,16,1,2,17,3,18,4,19,5,20,6:13,24,25)] svnames <- svnames[c(21,22,23,14,15,16,1,2,17,3,18,4,19,5,20,6:13,24,25)] list(qwstnum,svnames,svstuff) } ######################################################################################################### swaveqexPESTpdo <- function(dat) { # # Parameter estimation for seawaveqex model # clog <- dat[,"yobs"] centmp <- dat[,"icen"]=="c" # # The seasonal wave is computed using a period of 360 "days" # rather than 365 because it is required to have the period # be a multiple of 12. The following 3 lines define variables # that are used to "stretch" the period out to 365 days # tseas <- dat[,"sday"]/366 ipk <- floor(360*tseas) ipk[ipk==0] <- 1 # # Populate variables for julian day, mid-term and short-term # flow anomalies. and linear trend # jd <- dat[,"jday"] qamt <- dat[,"qanom30"] qast <- dat[,"qanom1"] lintnd <- dat[,"tnd"] # # Selection of the best seasonal wave model # The loop indeces correspond to the following cases: # - For wave class 1 (spmodcls=1), j is the choice of pulse input model # and k is the choice of the "half-life" (30, 60, 90, or 120 days) # - For wave class 2 (spmodcls=2) and k=1 or 2, j is the choice of pulse # input model from subclass A and k is the choice of "half-life" (30 or 60 days) # - For wave class 2 and k=3 or 4, j is the choice of pulse input model # from subclass B and k-2 is the choice of the "half-life" (30 or 60 days) #### find best wave from class 1 spmodcls <- 1 parxxx <- matrix(nrow=24,ncol=7) pvalxxx <- matrix(nrow=24,ncol=4) svcmaxt <- rep(NA,24) for (j in 1:6) { for (k in 1:4) { j2 <- (j-1)*4+k indcen <- !centmp # # For the specified wave model, the phase shift is selected from # 24 choices to maximize the likelihood function liksave <- rep(NA,24) for (imx in 1:24) { cmaxt <- imx/24 awave <- compwaveconvxx(cmaxt,j,k,spmodcls) wavest <- awave[ipk] tmpouta <- survreg(Surv(time=clog,time2=indcen,type="left") ~ wavest+qamt+qast+lintnd,dist="gaussian") liksave[imx] <- -tmpouta$loglik[2] } pckone <- order(liksave)[1] svcmaxt[j2] <- pckone/24 cmaxt <- pckone/24 awave <- compwaveconvxx(cmaxt,j,k,spmodcls) wavest <- awave[ipk] # # Now save the output for the wave model with the best phase shift # tmpouta <- survreg(Surv(time=clog,time2=indcen,type="left") ~ wavest+qamt+qast+lintnd,dist="gaussian") parxxx[j2,] <- c(summary(tmpouta)$table[1:5,1],tmpouta$scale,-tmpouta$loglik[2]) pvalxxx[j2,] <- c(summary(tmpouta)$table[2:5,4]) } } #### find best wave from class 2 spmodcls <- 2 parxxx2 <- matrix(nrow=24,ncol=7) pvalxxx2 <- matrix(nrow=24,ncol=4) svcmaxt2 <- rep(NA,24) for (j in 1:6) { for (k in 1:4) { j2 <- (j-1)*4+k indcen <- !centmp # # For the specified wave model, the phase shift is selected from # 24 choices to maximize the likelihood function liksave2 <- rep(NA,24) for (imx in 1:24) { cmaxt <- imx/24 awave <- compwaveconvxx(cmaxt,j,k,spmodcls) wavest <- awave[ipk] tmpouta <- survreg(Surv(time=clog,time2=indcen,type="left") ~ wavest+qamt+qast+lintnd,dist="gaussian") liksave2[imx] <- -tmpouta$loglik[2] } pckone2 <- order(liksave2)[1] svcmaxt2[j2] <- pckone2/24 cmaxt2 <- pckone2/24 awave <- compwaveconvxx(cmaxt2,j,k,spmodcls) wavest <- awave[ipk] # # Now save the output for the wave model with the best phase shift # tmpouta <- survreg(Surv(time=clog,time2=indcen,type="left") ~ wavest+qamt+qast+lintnd,dist="gaussian") parxxx2[j2,] <- c(summary(tmpouta)$table[1:5,1],tmpouta$scale,-tmpouta$loglik[2]) pvalxxx2[j2,] <- c(summary(tmpouta)$table[2:5,4]) } } # Select the best wave model for each class. If the regression coefficient for # the seasonal wave is negative do not choose that model because # the coefficient should be positive to make sense. Penalize class 2 liklihood # by adding 2. Compare likelihoods and choose best class. # likall <- (parxxx[,7]) ; likall[parxxx[,2]<0] <- NA pckone <- order(likall)[1] likval1 <- likall[pckone] likall2 <- (parxxx2[,7]) ; likall2[parxxx2[,2]<0] <- NA pckone2 <- order(likall2)[1] likval2 <- likall2[pckone2]+0 if (likval1 .005)) { sig0 <- sig1 neff <- 0 ssres <- 0 yimp <- y if(icen[1]) { yimp[1] <- sig0*mtmvnorm(mean=0,sigma=1,upper=y[1]/sig0)[[1]] } for (i in 2:ntot) { cxx <- rho^(jd[i]-jd[i-1])/del[i]/del[i-1] if (!icen[i]) { ztmp <- (y[i]-cxx*yimp[i-1])/(1-cxx^2)^0.5/sig0 if(ztmp>2.5) {ztmp <- 2.5 } else if (ztmp < (-2.5)) {ztmp <- (-2.5) } ssres <- ssres+ztmp^2*sig0^2 neff <- neff+1 yimp[i] <- ztmp*(1-cxx^2)^0.5*sig0 + cxx*yimp[i-1] } else { ztmp <- (y[i]-cxx*yimp[i-1])/(1-cxx^2)^0.5/sig0 if(ztmp<(-2.5)) {ztmp <- (-2.5) } neff <- neff-sig0^2*dnorm(ztmp)/pnorm(ztmp)*ztmp yimp[i] <- mtmvnorm(mean=0,sigma=1,upper=ztmp)[[1]] yimp[i] <- yimp[i]*(1-cxx^2)^0.5*sig0 + cxx*yimp[i-1] } } sig1 <- (ssres/neff)^0.5 # if (is.na(sig1)) {sig1 <- sig0} ###### if (is.na(sig1)) {sig1 <- sig0 } else if (sig1>1.01*sig0) {sig1 <- 1.01*sig0} } } sig1 } ######################################################################################################### evalmodlikxx <- function(parin,nrat,sig0,ytil,jd,icen,wavest,ifixs) { #### if ifixs=0, find updated process variance (signew) given rho, nrat, and starting value (sig0) rho <- exp(parin[2]) alph <- parin[1] if (ifixs==0) { signew <- estsigxx(rho,nrat,sig0,ytil,jd,icen,wavest,alph) } else { signew <- sig0 } #### compute conditional mimimum of -2logL omg <- (1+alph*(wavest)) del <- (1+nrat/omg^2)^0.5 y <- ytil/omg/del ntot <- length(y) sig0 <- signew cumlik <- 0 yimp <- y if(icen[1]) { yimp[1] <- sig0*mtmvnorm(mean=0,sigma=1,upper=y[1]/sig0)[[1]] } for (i in 2:ntot) { cxx <- rho^(jd[i]-jd[i-1])/del[i]/del[i-1] if (!icen[i]) { ztmp <- (y[i]-cxx*yimp[i-1])/(1-cxx^2)^0.5/sig0 if(ztmp>2.5) {ztmp <- 2.5 } else if (ztmp < (-2.5)) {ztmp <- (-2.5) } cumlik <- cumlik + ztmp^2 + 2*log(sig0) + log(1-cxx^2) yimp[i] <- ztmp*(1-cxx^2)^0.5*sig0 + cxx*yimp[i-1] } else { ztmp <- (y[i]-cxx*yimp[i-1])/(1-cxx^2)^0.5/sig0 if(ztmp<(-2.5)) {ztmp <- (-2.5) } cumlik <- cumlik-2*log(pnorm(ztmp)) yimp[i] <- mtmvnorm(mean=0,sigma=1,upper=ztmp)[[1]] yimp[i] <- yimp[i]*(1-cxx^2)^0.5*sig0 + cxx*yimp[i-1] } } cumlik <- cumlik + 2*sum(log(del[!icen]))+ 2*sum(log(omg[!icen])) + 500*(log(0.9) - log(rho))^4 + 0.5*(alph/2)^4 cumlik } ######################################################################################################### compwaveconvxx <- function(cmax,jmod,hlife,spmodcls) { # computes seasonal wave with cmax=time of peak, jmod=model choice (1:6), # and hlife=halflife, in months (1:4) # spmodcls=1 for "one-hump" models, spmodcls=2 for "two-hump" models del <- 1/12 txx <- seq(0,1,1/360) sxx <- txx if (spmodcls==1) { phi <- 12/hlife if(jmod==1) {wtx <- c(1,0,0,0,0,0,0,0,0,0,0,0); pkt <- 1/12} if(jmod==2) {wtx <- c(1,1,0,0,0,0,0,0,0,0,0,0); pkt <- 2/12} if(jmod==3) {wtx <- c(1,1,1,0,0,0,0,0,0,0,0,0); pkt <- 3/12} if(jmod==4) {wtx <- c(1,1,1,1,0,0,0,0,0,0,0,0); pkt <- 4/12} if(jmod==5) {wtx <- c(1,1,1,1,1,1,0,0,0,0,0,0); pkt <- 6/12} if(jmod==6) {wtx <- c(1,1,1,1,1,1,1,1,1,0,0,0); pkt <- 9/12} } else if (spmodcls==2) { if (hlife<=2) { phi <- 12/hlife if(jmod==1) {wtx <- c(1,1,0,0,1,1,1,1,0,0,0,0); pkt <- 8/12} if(jmod==4) {wtx <- c(1,1,1,0,0,1,1,1,0,0,0,0); pkt <- 8/12} if(jmod==3) {wtx <- c(1,1,1,1,0,0,1,1,0,0,0,0); pkt <- 4/12} if(jmod==2) {wtx <- c(1,1,0,0,0,1,1,1,1,0,0,0); pkt <- 9/12} if(jmod==5) {wtx <- c(1,1,1,0,0,0,1,1,1,0,0,0); pkt <- 9/12} if(jmod==6) {wtx <- c(1,1,1,1,0,0,1,1,1,1,0,0); pkt <- 10/12} } else { phi <- 12/(hlife-2) if(jmod==1) {wtx <- c(2,2,0,0,1,1,1,1,0,0,0,0); pkt <- 2/12} if(jmod==4) {wtx <- c(1.5,1.5,1.5,0,0,1,1,1,0,0,0,0); pkt <- 3/12} if(jmod==3) {wtx <- c(1,1,1,1,0,0,2,2,0,0,0,0); pkt <- 8/12} if(jmod==2) {wtx <- c(2,2,0,0,0,1,1,1,1,0,0,0); pkt <- 2/12} if(jmod==5) {wtx <- c(1.5,1.5,1.5,0,0,0,1,1,1,0,0,0); pkt <- 3/12} if(jmod==6) {wtx <- c(1.5,1.5,1.5,1.5,0,0,1,1,1,1,0,0); pkt <- 4/12} } } rho <- exp(-phi) z0xx <- rho^(txx) r12 <- rho^(-del*c(1:12)) con <- wtx[1]*(r12[1]-1) for(k in 2:12) {con <- con+wtx[k]*(r12[k]-r12[k-1])} con <- rho/(1-rho)*con z0xx <- z0xx*con zmat <- matrix(nrow=length(txx),ncol=12) pckm <- matrix(nrow=length(txx),ncol=12) ntot <- length(txx) pckm[,1] <- c(txx<=del) for (k in 2:12) {pckm[,k] <- c(txx>(k-1)*del & txx<=k*del)} zmat[,1] <- rho^txx*(rho^(-replace(txx,txx>1/12,1/12))-1) for (k in 2:12) { ztmp <- rep(0,ntot) for (j in 1:(k-1)) {ztmp[pckm[,j]] <- 0} ztmp[pckm[,k]] <- 1-rho^(txx[pckm[,k]]-(k-1)*del) if(k<12){ for (j in (k+1):12) {ztmp[pckm[,j]] <- rho^(txx[pckm[,j]]-del*k)-rho^(txx[pckm[,j]]-(k-1)*del)} } zmat[,k] <- ztmp } sst <- z0xx for (k in 1:12) { sst <- sst+wtx[k]*zmat[,k]} sst <- sst/phi medxx <- (max(sst)+min(sst))/2 rngxx <- 2*(max(sst)-medxx) sst <- (sst-medxx)/rngxx ###################### if(cmax<=pkt) { txx2 <- (txx+1-pkt+cmax)} if(cmax>pkt) {txx2 <- txx-pkt+cmax} txx2[txx2>1] <- txx2[txx2>1]-1 otmp <- order(txx2); sst <- sst[otmp] sst } ######################################################################################################### swaveqexCSIM <- function(ngen,parxxx,datobs,datall,spmodcls) { # generates conditional simulations from seawaveqex model ## define seasonal wave parameters etc hlife <- parxxx[10] ; cmax <- parxxx[8] ; modn <- parxxx[9] swave <- compwaveconvxx(cmax,modn,hlife,spmodcls) sdum <- (1:365)/366 ipk <- floor(360*sdum) ipk[ipk==0] <- 1 swave <- swave[ipk] waveobs <- swave[datobs[,"sday"]] waveall <- swave[datall[,"sday"]] datobs[,"swave"] <- waveobs datall[,"swave"] <- waveall ## define other model parameters etc. alph <- parxxx[12] ; rho <- parxxx[13] ; sig <- parxxx[14] omgobs <- (1+alph*(0.0+waveobs)) omgall <- (1+alph*(0.0+waveall)) int0 <- parxxx[1] ; wc0 <- parxxx[2]; mtc0 <- parxxx[3]; stc0 <- parxxx[4]; tndc0 <- parxxx[5] muobs <- int0*datobs[,"int"]+wc0*datobs[,"swave"]+mtc0*datobs[,"qanom30"]+stc0*datobs[,"qanom1"]+tndc0*datobs[,"tnd"] muall <- int0*datall[,"int"]+wc0*datall[,"swave"]+mtc0*datall[,"qanom30"]+stc0*datall[,"qanom1"]+tndc0*datall[,"tnd"] yhatobs <- mtc0*datobs[,"qanom30"]+stc0*datobs[,"qanom1"]+tndc0*datobs[,"tnd"] yhatall <- mtc0*datall[,"qanom30"]+stc0*datall[,"qanom1"]+tndc0*datall[,"tnd"] yobs <- datobs[,"yobs"] ; ytilobs <- (yobs-muobs)/omgobs yall <- datall[,"yobs"] ; ytilall <- (yall-muall)/omgall jdayobs <- datobs[,"jday"] jdayall <- datall[,"jday"] icen <- datobs[,"icen"]=="c" svcondsim <- matrix(nrow=length(jdayall),ncol=ngen) for (igen in 1:ngen) { # impute censored values if necessary yimp <- ytilobs if (sum(icen) > 0) { yimp <- impcenvals(rho,sig,yimp,icen,jdayobs) } # now do conditional simulation ysim <- condsim(rho,sig,yimp,jdayobs,jdayall) ysim <- ysim*omgall + muall yimp <- yimp*omgobs + muobs svcondsim[,igen] <- ysim } ## compute various statistics etc for output ... svcondmean <- apply(10^svcondsim,1,mean) svcondmean <- log10(svcondmean) svcondsim1 <- svcondsim[,1] nall <- length(jdayall) nyrs <- round(nall/365,0) tmpannmax <- matrix(nrow=nyrs,ncol=ngen) for (i in 1:ngen) { tmpyy <- 10^svcondsim[,i] tmpannmax[,i] <- apply(matrix(tmpyy,nrow=365),2,max) } maxout <- apply(tmpannmax,1,mean) modoutobs <- data.frame(yhatobs,muobs,omgobs,svcondsim1[jdayobs]) names(modoutobs) <- c("yhat","mu","omg","condsim") modoutall <- data.frame(yhatall,muall,omgall,svcondsim1,svcondmean) names(modoutall) <- c("yhat","mu","omg","condsim","condmean") list(maxout,modoutobs,modoutall,svcondsim) } ######################################################################################################### impcenvals <- function(rho,sig,yimp,icen,jdayobs) { # # generates conditional values for censored data only # nobs <- length(yimp) ncum <- 1:nobs indu <- !icen ; iu <- ncum[indu] nu <- sum(indu) nblkc <- nu+1 ncb <- rep(0,nblkc) if(icen[1]) {ncb[1] <- iu[1]-1} for (i in 2:nu) {ncb[i] <- iu[i]-iu[i-1]-1} if (icen[nobs]) {ncb[nu+1] <- nobs-iu[nu]} t1 <- jdayobs #### ready to impute new values for censored data for (i in 1:(nu+1)) { if(i==1 & ncb[1]>0) { yblk <- yimp[1:ncb[1]] yends <- yimp[ncb[1]+1] mtmp <- matrix(nrow=ncb[1]+1,ncol=ncb[1]+1) for (m in 1:(ncb[1]+1)) { for (k in m:(ncb[1]+1)) { if(k==m) {mtmp[k,k] <- 1} if(k>m) { mtmp[m,k] <- rho^(t1[k]-t1[m]) mtmp[k,m] <- mtmp[m,k] } } } mtmpx <- as.matrix(mtmp[1:ncb[1],1:ncb[1]]); gtmpx <- as.matrix(mtmp[1:ncb[1],ncb[1]+1],ncol=1) mtmpy <- mtmp[ncb[1]+1,ncb[1]+1] cvar <- mtmpx-gtmpx%*%t(gtmpx)/mtmpy cvarinv <- solve(cvar) cmu <- gtmpx%*%as.matrix(yends,ncol=1)/mtmpy yimp[1:ncb[1]] <- rtmvnorm(n=1,mean=c(cmu),,sigma=sig^2*cvar,upper=yblk,algorithm='gibbs',burn.in.samples=100) } if (i>1 & i<(nu+1) & ncb[i]>0) { yblk <- yimp[(iu[i-1]+1):(iu[i]-1)] yends <- yimp[c(iu[i-1],iu[i])] mtmp <- matrix(nrow=ncb[i]+2,ncol=ncb[i]+2) for (m in 1:(ncb[i])) { for (k in m:(ncb[i])) { if(k==m) {mtmp[k,k] <- 1} if(k>m) { mtmp[m,k] <- rho^(t1[iu[i-1]+k]-t1[iu[i-1]+m]) mtmp[k,m] <- mtmp[m,k] } } } for (m in 1:ncb[i]) { mtmp[m,ncb[i]+1] <- rho^(t1[iu[i-1]+m]-t1[iu[i-1]]) mtmp[ncb[i]+1,m] <- mtmp[m,ncb[i]+1] mtmp[m,ncb[i]+2] <- rho^(t1[iu[i]]-t1[iu[i-1]+m]) mtmp[ncb[i]+2,m] <- mtmp[m,ncb[i]+2] } mtmp[ncb[i]+1,ncb[i]+1] <- 1 mtmp[ncb[i]+2,ncb[i]+2] <- 1 mtmp[ncb[i]+1,ncb[i]+2] <- rho^(t1[iu[i]]-t1[iu[i-1]]) mtmp[ncb[i]+2,ncb[i]+1] <- mtmp[ncb[i]+1,ncb[i]+2] mtmpx <- as.matrix(mtmp[1:ncb[i],1:ncb[i]]); gtmpx <- matrix(mtmp[1:ncb[i],c(ncb[i]+1,ncb[i]+2)],ncol=2) mtmpy <- mtmp[c(ncb[i]+1,ncb[i]+2),c(ncb[i]+1,ncb[i]+2)] cvar <- mtmpx-gtmpx%*%solve(mtmpy)%*%t(gtmpx) cmu <- gtmpx%*%solve(mtmpy)%*%as.matrix(yends,ncol=1) cvarinv <- solve(cvar) yimp[(iu[i-1]+1):(iu[i-1]+ncb[i])] <- rtmvnorm(n=1,mean=c(cmu),sigma=sig^2*cvar,upper=yblk,algorithm='gibbs',burn.in.samples=100) } if(i==(nu+1) & ncb[nu+1]>0) { yblk <- yimp[(iu[nu]+1):(iu[nu]+ncb[nu+1])] yends <- yimp[iu[nu]] mtmp <- matrix(nrow=ncb[nu+1]+1,ncol=ncb[nu+1]+1) for (m in 1:(ncb[nu+1]+1)) { for (k in m:(ncb[nu+1]+1)) { if(k==m) {mtmp[k,k] <- 1} if(k>m & k<(ncb[nu+1]+1)) { mtmp[m,k] <- rho^(t1[iu[nu]+k]-t1[iu[nu]+m]) mtmp[k,m] <- mtmp[m,k] } if(k>m & k==(ncb[nu+1]+1)) { mtmp[m,k] <- rho^(t1[iu[nu]+m]-t1[iu[nu]]) mtmp[k,m] <- mtmp[m,k] } } } mtmpx <- as.matrix(mtmp[1:ncb[nu+1],1:ncb[nu+1]]); gtmpx <- as.matrix(mtmp[1:ncb[nu+1],ncb[nu+1]+1],ncol=1) mtmpy <- mtmp[ncb[nu+1]+1,ncb[nu+1]+1] cvar <- mtmpx-gtmpx%*%t(gtmpx)/mtmpy cmu <- gtmpx%*%as.matrix(yends,nrow=1)/mtmpy cvarinv <- solve(cvar) yimp[(iu[nu]+1):(iu[nu]+ncb[nu+1])] <- rtmvnorm(n=1,mean=c(cmu),sigma=sig^2*cvar,upper=yblk,algorithm='gibbs',burn.in.samples=100) } } yimp } ######################################################################################################### condsim <- function(rho,sig,yimp,jdayobs,jdayall) { # # generates conditional simulations # nobs <- length(yimp) nall <- length(jdayall) ysim <- rep(NA,nall) tobs <- jdayobs tall <- jdayall for (i in 1:nobs) { ysim[tobs[i]] <- yimp[i] if(i==1 & tobs[1]>1) { n <- tobs[1] nm1 <- n-1 yends <- yimp[1] mtmp <- matrix(nrow=n,ncol=n) for (m in 1:n) { for (k in m:n) { if(k==m) {mtmp[k,k] <- 1} if(k>m) { mtmp[m,k] <- rho^(k-m) mtmp[k,m] <- mtmp[m,k] } } } mtmpx <- as.matrix(mtmp[1:nm1,1:nm1]); gtmpx <- as.matrix(mtmp[1:nm1,n]) mtmpy <- mtmp[n,n] cvar <- mtmpx-gtmpx%*%t(gtmpx)/mtmpy cvarinv <- solve(cvar) cmu <- gtmpx%*%as.matrix(yends,ncol=1)/mtmpy ysim[1:nm1] <- rmvnorm(n=1,mean=cmu,sigma=sig^2*cvar,method="chol") } if (i>1 & i<=nobs & (tobs[i]-tobs[max(1,i-1)])>0) { n <- tobs[i]-tobs[i-1]+1 nm2 <- n-2 yends <- yimp[c(i-1,i)] mtmp <- matrix(nrow=n,ncol=n) for (m in 1:nm2) { for (k in m:nm2) { if(k==m) {mtmp[k,k] <- 1} if(k>m) { mtmp[m,k] <- rho^(k-m) mtmp[k,m] <- mtmp[m,k] } } } for (m in 1:nm2) { mtmp[m,nm2+1] <- rho^(m) mtmp[nm2+1,m] <- mtmp[m,nm2+1] mtmp[m,nm2+2] <- rho^(nm2-m+1) mtmp[nm2+2,m] <- mtmp[m,nm2+2] } mtmp[nm2+1,nm2+1] <- 1 mtmp[nm2+2,nm2+2] <- 1 mtmp[nm2+1,nm2+2] <- rho^(tobs[i]-tobs[i-1]) mtmp[nm2+2,nm2+1] <- mtmp[nm2+1,nm2+2] mtmpx <- as.matrix(mtmp[1:nm2,1:nm2]); gtmpx <- matrix(mtmp[1:nm2,c(nm2+1,nm2+2)],ncol=2) mtmpy <- mtmp[c(nm2+1,nm2+2),c(nm2+1,nm2+2)] cvar <- mtmpx-gtmpx%*%solve(mtmpy)%*%t(gtmpx) cmu <- gtmpx%*%solve(mtmpy)%*%as.matrix(yends,ncol=1) ysim[(tobs[i-1]+1):(tobs[i]-1)] <- rmvnorm(n=1,mean=c(cmu),sigma=sig^2*cvar,method="chol") } if(i==nobs & tobs[i]m) { mtmp[m,k] <- rho^(k-m) mtmp[k,m] <- mtmp[m,k] } } } for (m in 1:nm1) { mtmp[m,n] <- rho^(m) mtmp[n,m] <- mtmp[m,n] } mtmp[n,n] <- 1 mtmpx <- as.matrix(mtmp[1:nm1,1:nm1]); gtmpx <- as.matrix(mtmp[1:nm1,n]) mtmpy <- mtmp[n,n] cvar <- mtmpx-gtmpx%*%t(gtmpx)/mtmpy cvarinv <- solve(cvar) cmu <- gtmpx%*%as.matrix(yends,ncol=1)/mtmpy ysim[(tobs[i]+1):nall] <- rmvnorm(n=1,mean=cmu,sigma=sig^2*cvar,method="chol") } } ysim } ######################################################################## that's it!