Nothing
#library(devtools)
#library(roxygen2)
# Nchars1 = unlist(regexec("x",getwd())) + 1
# file0 = paste(substr(getwd(),start = 1, stop = Nchars1),"myprofile.r", sep = "" )
# FilesToLoad = c(file0)
# sapply(FilesToLoad,source,.GlobalEnv)
#------------ package main -----------------
#setwd('C:/Users/MP-User/Dropbox/E2/Webfolder/posts')
# coll <- c(wes.palette(4, "Moonrise2")[4],wes.palette(4, "Moonrise2")[1:3],
# wes.palette(4, "GrandBudapest")[1:4])
# @keywords K
# @seealso
# @aliases
coll <- c("#29211F","#798E87", "#C27D38", "#CCC591", "#F1BB7B","#FD6467",
"#5B1A18","#D67236") # fancy colors from wesanderson package
lwd1 <- fontt <- 2
cexx=1.3
# parr <- par(mfrow = c(2,2), bg=wes.palette(4, "Chevalier")[3], bty ="n", fg = 'black' ,
# font.lab = fontt,
# font.axis = fontt, font.main = fontt, col.axis = 'black',col.lab = 'black',
# cex.axis = 1.4, pch = 21, tck = -0.02, cex.lab = 1.6 ,cex.main = 2,las=1,
# mar = c(5, 4, 4, 2) + 0.4)
#' plott
#'
#' Plotting longitudinal series.
#'
#' Sets default parameters to get a nicer figure. If y is given then a scatter
#' plot is created. y can also be of class Date. If add=TRUE, add the series to
#' an existing device. Make a sensible choice as to which series should be
#' plotted first so that the ylim and xlim cover are sufficient.
#'
#' @param x series to be plotted.
#' @param y possible second series, if provided a scatter plot is created.
#' @param add if \code{add = TRUE} the series is added to existing active device. The active device graphical
#' parameters must match, meaning it must be created using \code{return.to.default=F}.
#' @param pch same as in \code{\link{par}}.
#' @param xlab a label for the x axis, same as in \code{\link{plot.default}}.
#' @param main main title of the plot, same as in \code{\link{plot.default}}.
#' @param col Color of the series, same as in \code{\link{par}}.
#' @param ty character indicating the type of plotting, any of the types as in \code{\link{plot.default}}.
#' @param return.to.default for reverting back to previous par settings.
#' Default is \code{return.to.default=TRUE}.
#' @param ... more graphical parameters can be given as arguments.
#' @return Called for its side effect.
#' @seealso \code{\link{par}}, \code{\link{plot.default}}
#' @export plott
plott = function(x,y=c(1:length(x)),add=FALSE,pch=19,xlab="",col=1,
main=NULL,ty="b",return.to.default=TRUE,...){
lwd1 <- fontt <- 2
if(add){
lines(x~y,xlab="",ylab = "",ty=ty,lwd=lwd1,col=col,pch=pch)
}
if(!add){
pardefault <- par(no.readonly = TRUE)
parr = par(bg="white", bty ="n", fg = 'black' , font.lab = fontt,
font.axis = fontt, font.main = fontt, col.axis = 'black',col.lab = 'black',
cex.axis = 1.4, pch = 21, tck = -0.02, cex.lab = 1.6 ,cex.main = 2,las=1,
mar = c(5, 4, 4, 2) + 0.4)
plot(x~y,xlab=xlab, ylab = "",ty=ty,lwd=2,main=main,pch=pch,col=col,...)
grid(col = "grey")
if(return.to.default) { par(pardefault)}
} }
#' tsideplot
#'
#' Create a plot of two series vertical axes on both left and right side.
#'
#' Create a plot of two series with y-axes on both left and right side.
#' \code{Set return.to.default=TRUE} to keep the new settings, otherwise
#' default to revert to previous par values. xaxis parameter is the optional
#' xaxis, if not provided then \code{if(is.null(xaxis)) {xaxis=
#' c(1:length(series1))}} is used.
#'
#' @param series1,series2 First and second series to be plotted.
#' @param main main title of the plot, same as in \code{\link{plot.default}}.
#' @param return.to.default for reverting back to previous par settings.
#' @param xaxis Optional, the xaxis to be used, see details.
#' @param col Color of the second series, same as in \code{\link{par}}.
#' @param ... more graphical parameters can be given as arguments.
#' @return Called for its side effect.
#' @seealso \code{\link{par}}, \code{\link{plot.default}}.
#' @export tsideplot
tsideplot <- function(series1,series2,main="", return.to.default=T,
xaxis=NULL, col='red',...){
if (NROW(series1) != NROW(series2) ) { stop("Problem is that NROW(series1) != NROW(series2)") }
lwd1 <- fontt <- 2
pardefault <- par(no.readonly = TRUE)
par(mfrow = c(1,1),bg="white", bty ="n", fg = 'black' , font.lab = fontt,
font.axis = fontt, font.main = fontt, col.axis = 'black',col.lab = 'black',
cex.axis = 1.4, pch = 21, tck = -0.02, cex.lab = 1.6 ,cex.main = 2,las=1,
mar = c(5, 4, 4, 4) + 0.1)
if(is.null(xaxis)) {xaxis= c(1:length(series1))}
plot(series1~xaxis,ty="b",pch=19,xlab="",ylab="",col = 1,lwd=lwd1,...)
par(new=TRUE)
plot(series2~xaxis,main=main,bty="o",ty="b",pch=19,xlab="",ylab="",
axes=F,col = col,lwd=lwd1,...)
axis(4,at=pretty(series2),col=col,col.ticks=col,col.axis=col)
grid(col = "grey")
if(return.to.default){par(pardefault)}
}
#' mplott
#'
#' Multivariate plot.
#'
#' Multivariate plot. Limited to 5 series. Legend is added automatically using
#' \code{colnames(x)}.
#'
#' @param x a matrix to be plotted
#' @param wherelegend where to place the legend
#' @param textlegend what should the legend read (see details).
#' @param main main title of the plot, same as in \code{\link{plot.default}}.
#' @param return.to.default for reverting back to previous par settings.
#' Default is \code{return.to.default=TRUE}
#' @param ... more graphical parameters can be given as arguments.
#' @return called for its side effect.
#' @seealso \code{\link{par}}, \code{\link{plot.default}}.
#' @export mplott
mplott <- function(x,wherelegend='bottomleft',textlegend=colnames(x),main="",
return.to.default=T,...){
lwd1 <- fontt <- 2
pardefault <- par(no.readonly = TRUE)
l = NCOL(x)
pchh <- c(19,1:4)
if (l>5) {
l=5
textlegend=colnames(x)[1:5]
print("Only columns 1:5 are plotted")
}
par(bg="white", bty ="n", fg = 'black' , font.lab = fontt,
font.axis = fontt, font.main = fontt, col.axis = 'black',col.lab = 'black',
cex.axis = 1.4, pch = 21, tck = -0.02, cex.lab = 1.6 ,cex.main = 2,las=1,
mar = c(5, 4, 4, 4) + 0.1)
x <- as.matrix(x)
plot(x[,1],ty = "b",ylim = c(min(na.omit(x)),max(na.omit(x))),xlab = "",
ylab = "",lwd=2,pch=pchh[1],main=main,...)
for (i in 2:l){
if(l==1) {next}
lines(x[,i],col=i,ty = "b",lty=i,lwd=2,pch=pchh[i],...)
}
if(is.null(textlegend)) { textlegend = paste("Series", 1:l) }
legend(wherelegend,textlegend,col=1:l,lty=1:l,cex=1.3,text.font=2,
ncol=2,bty="n",text.col=1:l,lwd = rep(2,l),pch=pchh[1:l])
grid(col='grey')
if(return.to.default){par(pardefault)}
}
#mplott(x,wherelegend='topleft')
#x <- cbind(rnorm(100),rnorm(100),rnorm(100))
# mplott(cbind(x,y,z))
#' lagmat
#'
#' Creates a lagged matrix with the desired number of lags.
#'
#' @param x the series to be lagged
#' @param lags number of lags desired
#' @return matrix with dimension \code{[NROW(x),length(lags)]}
#' @examples
#'
#' x = rnorm(100)
#' lx <- lagmat(x,2)
#' tail(lx)
#' tail(x)
#'
#' @export lagmat
lagmat <- function(x,lags){
TT <- NROW(x) ; colnam <- NULL
lagmatt <- matrix(nrow = TT, ncol = lags)
for (i in 1:lags){
lagmatt[(i+1):(TT),i] <- x[1:(TT-i)]
colnam[i] <- paste('Lag',i,sep="")
}
lagmatt[1:lags,] <- NA
colnames(lagmatt) <- colnam
lagmatt
}
#' linpred
#'
#' Provides linear regression based predictions from a \code{y~x} type model
#' using recursive or rolling regression.
#'
#' The training is done using the direct method: \eqn{y_{1 : (t+h-1)} = \beta
#' x_{1:(t-1)} + \varepsilon_{1:(t+h-1)} } and the forecast is made at time
#' (t+h) as \eqn{\widehat{y}_{t+h} = \widehat{\beta} x_t}.
#'
#' @param y a series to be predicted
#' @param x a numeric or matrix of explanatory variables
#' @param h The horizon for which you would like to have the prediction for
#' (see details)
#' @param wind the size of the rolling window or the initial training period if
#' recursive is used
#' @param rr recursive or rolling window? Possible values are
#' \code{c("Rec","Rol")}
#' @return vector of prediction values with the same dimension as the original
#' series. The first \code{wind} values are NA's
#' @examples
#'
#' x = rnorm(100)
#' lx <- lagmat(x,2)
#' tail(lx)
#' tail(x)
#' out <- linpred(x,lx)
#' plott(x, return.to.default=FALSE)
#' plott(out,add=TRUE,col=2)
#'
#' @export linpred
linpred <- function(y,x,h=1,wind=NULL,rr=c("Rec")){
if (rr != "Rol" & rr !="Rec") { stop("rr is not properly defined") }
if (NROW(y) != NROW(x) ) { stop("Problem is that NROW(y) != NROW(x)") }
if(is.null(wind)){wind=0.25*length(y)}
predd <- NULL
TT <- length(y)
x=as.matrix(x) # if x is a numeric make it a matrix
if(rr=='Rec'){
for (j in (1+wind):TT){
lm0 <- lm(y[1:(j-h)]~x[1:(j-h),]) # train the coefficietns, stop before i+wind-h.
predd[j] <- lm0$coef%*%c(1,x[(j-h),])
}}
if(rr=='Rol'){
for (j in 1:(TT-wind)){
lm0 <- lm(y[j:(wind+j-h)]~x[j:(wind+j-h),]) # train the coefficietns, stop before i+wind-h.
predd[(j+wind)] <- lm0$coef%*%c(1,x[(j+wind-h),])
} }
predd
}
#' FCIplot
#'
#' Estimate and plot prediction standard deviation. Given the series, the
#' function estimate point prediction based on AR(1) model and, using the
#' resdiuals from this simple model, estimate an ARCH model to estimate the
#' prediction standard deviation. If \code{plott=TRUE}, a plot of the most
#' recent \code{k} values is created.
#'
#' Estimate and plot prediction confidence intervals based on AR-ARCH model.
#'
#' @param series series to be plotted.
#' @param plott should a plot be created? default is \code{plott=TRUE}.
#' @param wind1 window size for the AR component (see details).
#' @param wind2 window size for the ARCH component (see details).
#' @param k if \code{plott=TRUE}, \code{tail(series,k)} will be plotted.
#' @param rrr1 will the AR model be estimated using Recursive ("Rec") or
#' Rolling ("Rol") window?
#' @param rrr2 will the ARCH model be estimated using Recursive ("Rec") or
#' Rolling ("Rol") window?
#' @param main main title of the plot, same as in \code{\link{plot.default}}.
#' @return vector of prediction's standard deviation.
#' @examples
#'
#' par(mfrow = c(2,1))
#' out <- FCIplot(rnorm(100),plott=TRUE,k=30)
#' plott(out,main="The out-of-sample standard deviation")
#'
#' @export FCIplot
FCIplot <- function(series,plott=TRUE,wind1=24,wind2=60,k=60,rrr1="Rec",rrr2="Rec",main="series"){
if (wind2< wind1) { stop("wind2 has to be > than wind1") }
TT <- length(series)
temp <- linpred(series[2:TT],series[1:(TT-1)],h=1, wind=wind1,rr=rrr1)
ertemp <- (series[2:TT]-temp)^2
TT2 <- length(ertemp)
sdestimate <- linpred(ertemp[2:TT2], ertemp[1:(TT2-1)],
wind = wind2, rr = rrr2)
# We lost 1 observation to construct the AR forecast and 1 to construct the ARCH
up1 = series[(wind2+2):(TT-1)] + 1*sqrt(abs(na.omit(sdestimate)))
down1 = series[(wind2+2):(TT-1)] - 1*sqrt(abs(na.omit(sdestimate)))
up2 = series[(wind2+2):(TT-1)] + 2*sqrt(abs(na.omit(sdestimate)))
down2 = series[(wind2+2):(TT-1)] - 2*sqrt(abs(na.omit(sdestimate)))
if(k>length(down2)){k=length(down2) ; warning( "k was set to a lower value" )}
if(plott){
plot(tail(series,k),ylim=c(min(tail(down2,k)),max(tail(up2,k))),ty="b",pch=19,xlab="",
main=main,ylab="")
grid()
polygon(x = c(rev(1:k),(1:k)), y = c(rev(tail(down2,k)),tail(up2,k)),
col = rgb(0,0,0,0.2),border ='red')
polygon(x = c(rev(1:k),(1:k)), y = c(rev(tail(down1,k)),tail(up1,k)),
col = rgb(0,0,0,0.2),border ='blue')
}
return(sdestimate^.5)
}
Any scripts or data that you put into this service are public.
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.