############################################################### ######### HALF-SIBLING AUTOREGRESSION : A SIMULATION ########## ############################################################### ############################################### ######### Loading Required Packages ########### ############################################### require(RandomFields) require(glmnet) require(zoo) ############################################### ###### General Setting and Parameters ######### ############################################### lPeriod = 200 # length of 1 period ltransit = .05*lPeriod # length of 1 transit SingleTransit = rep(c(0,-.25),times = c(lPeriod-ltransit,ltransit)) # F_t over 1 period alpha = 2*ltransit+1 delta = .03*lPeriod ############################################### ################# Training #################### ############################################### nPeriods = 10 # Number of training periods nGen = lPeriod * nPeriods + 2*(alpha + delta) # number of elements to generate, such that after lagging the time series, we get exactly nPeriods periods. scaleCov = .4 # Scale of covariance for the Trend,expressed with unit = lPeriod (= length of 1 period) model <- RMgneiting(var = 1.3 , scale = scaleCov) Trend = RFsimulate(model, x=seq(0, nGen/lPeriod,length=nGen)) Noise = rnorm(nGen) Transit = zoo(rep(SingleTransit,length.out = nGen)) ytrain = zoo(.05*as.vector(Trend) + .01*Noise + Transit) # zoo makes lagging of time series easy... plot(1:length(ytrain)/lPeriod, ytrain, type='l', xlab = 'time') xytrain = data.frame(sapply(c(0,(-delta):(-1)-alpha, 1:delta + alpha), function(i) as.numeric(lag(ytrain,i,na.pad=TRUE)))) xytrain = na.omit(xytrain) # lagging introduces NA values at the edges. Here we suppress them. ytrain = xytrain[,1] xtrain = I(as.matrix(xytrain[,-1])) LS = glmnet(xtrain,ytrain, alpha = 0, lambda = 0.1, intercept = TRUE) ############################################### ################# Testing ##################### ############################################### nTestPeriods = 4 TrendTest = RFsimulate(model, x=seq(0, nTestPeriods,length=nTestPeriods*lPeriod)) plot(TrendTest, type='l') NoiseTest = rnorm(nTestPeriods*lPeriod) TransitTest = rep(SingleTransit,nTestPeriods) ytest = zoo(.05*as.vector(TrendTest) + .01*NoiseTest + TransitTest) dev.off(); plot(1:length(ytest)/lPeriod, ytest, type='l', xlab = 'time') xytest = data.frame(sapply(c(0,(-delta):(-1)-alpha, 1:delta + alpha), function(i) as.numeric(lag(ytest,i,na.pad=TRUE)))) xytest = na.omit(xytest) ytest = xytest[,1] xtest = I(as.matrix(xytest[,-1])) yhat_test = predict(LS, newx = xtest) # yhat_test = rowMeans(xtest) # Nice Comparison ############################################### ############### Published Plot ################ ############################################### par(mar=c(2,1,1,1)+0.1) # margin for publication plot(1:length(ytest)/lPeriod, -.5 + TransitTest[(delta+alpha) + 1:length(ytest)] ,xlab = '', ylab='',yaxt = 'n', ylim = c(-.8, 1.25), type='l') lines(1:length(ytest)/lPeriod,ytest) lines(1:length(ytest)/lPeriod,.5+yhat_test, col = 'blue') lines(1:length(ytest)/lPeriod,1+(ytest - yhat_test), col='blue') # ############################################### # ######## Other (unpublished) Plots ############ # ############################################### # # plotIndices = 1:(min(length(ytrain),5*lPeriod)); plotTime = plotIndices/lPeriod # Plotting Range # # # Stationary Covariance of the Trend # dev.off(); plot(model, xlim = c(0,2.5*scaleCov)) # # # Plot of the Trend # dev.off(); plot(Trend) # # # First 5 Periods of the Training Data # dev.off(); plot(plotTime, ytrain[plotIndices],type='l', xlab = 'time', ylab = 'ytrain') # # # Residuals on Training Set # dev.off(); plot(plotTime , predict(LS,xtrain[plotIndices,]) - ytrain[plotIndices], type='l', xlab = 'time', ylab= 'residuals', main = 'Residuals on training set') # # # Ridge Regression coefficients # dev.off(); matplot(coef(LS),type = 'l', lty = 'solid') # # # ytrain and regressors xtrain = ytrain_{t+W} # matplot(plotTime,xytrain[plotIndices,], type='l', lty='solid')