utils::globalVariables(c("aaply"))
logLikeFn<-function(residual){
SS =sum(residual^2, na.rm=TRUE)
n =sum(!is.na(residual))
se =(SS/n)^.5
#res =(log(1/(2*pi))-n*log(se)-SS/(2*se^2))/2
res=n/2*log(3.14159265359*2)+n*log(se)+SS/(2*se*se)
return(res)}
objLog=function(residual){
dmn=list(params=c("ss","lav","ll","n"),
iter =dimnames(residual)$iter)
rtn=FLPar(array(as.numeric(NA),dim=unlist(laply(dmn, length)),dimnames=dmn))
rtn[ "ss",]=aaply(residual%*%residual,6,sum,na.rm=TRUE)
rtn["lav",]=aaply(abs(residual), 6,sum,na.rm=TRUE)
rtn[ "ll",]=aaply(residual, 6,function(x) logLikeFn(x))
rtn[ "n",]=aaply(residual, 6,function(x) sum(!is.na(x)))
rtn}
calcQFLQuant<-function(stock,index,when=0.5,na.rm=T){
when =pmin(pmax(when,0),1)
stock=stock[,-dim(stock)[2]]*(1-when)+stock[,-1]*when
yrs=dimnames(index)$year
yrs=yrs[yrs%in%dimnames(stock)$year]
index=index[,yrs]
stock=stock[,yrs]
q =exp(apply(log(index)%-%log(stock),6,sum, na.rm=na.rm)/dim(index)[2])
#sigma=calcSigma(log(index),log(stock%*%q))
return(q)}
calcObjFn=function(object,index="missing",when=0.5,calcq=TRUE,na=10e6){
if (missingArg(index))
index=FLQuants(index(object,FALSE))
else if (is.FLQuant(index))
index=FLQuants(index)
if (!calcq)
q=params(object)[grep("q",dimnames(params(object))$params)]
else
q=FLPar(q=laply(index,function(x) calcQFLQuant(stock(object),x)))
res=FLPars(mlply(data.frame(component=seq(length(index))), function(component){
obs=index[[component]]%/%q[component]
hat=stock(object,when=when)
yrs=dimnames(obs)$year
yrs=yrs[yrs%in%dimnames(hat)$year]
residuals=log(obs[,yrs]/hat[,yrs])
residuals[is.na(residuals)]=na
objLog(residuals)}))
attributes(res)=NULL
res}
calcObjFn2=function(object,when=0.5,calcq=TRUE,na=10e6){
if (missingArg(index))
index=FLQuants(index(object,FALSE))
else if (is.FLQuant(index))
index=FLQuants(index)
if (!calcq)
q=params(object)[grep("q",dimnames(params(object))$params)]
else
q=FLPar(q=laply(index,function(x) calcQFLQuant(stock(object),x)))
res=FLPars(mlply(data.frame(component=seq(length(index))), function(component){
obs=index[[component]]%/%q[component]
hat=stock(object,when=when)
yrs=dimnames(obs)$year
yrs=yrs[yrs%in%dimnames(hat)$year]
residuals=log(obs[,yrs]/hat[,yrs])
residuals[is.na(residuals)]=na
objLog(residuals)}))
attributes(res)=NULL
res}
Add the following code to your website.
For more information on customizing the embed code, read Embedding Snippets.