Nothing
# nscosinor.initial.R
# Estimate starting values for seasonal variance based on a linear trend and given values for tau
# Dec 2011
# set up variables and matrices
nscosinor.initial=function(data,response,tau,lambda=1/12,n.season){
attach(data, warn.conflicts = FALSE)
on.exit(detach(data))
k<-1; # Assume just one season
kk<-2*(k+1);
n<-nrow(data);
F<-rep(c(1,0),k+1);
alpha_j<-matrix(0,kk,n+1);
G<-matrix(0,kk,kk);
G[1,1]=1; G[1,2]=lambda; G[2,2]=1;
# linear model
time=1:n
model=glm(response~time,data=data)
## put predictions into alpha
# 1. trend
alpha_j[1,2:(n+1)]=fitted(model)
# 2. season
sdr=sd(resid(model))
alpha_j[3,2:(n+1)]=(sdr/10)*cos(2*pi*(1:n+1)/12) # sinusoid with amplitude equal to 10% of standard deviation of residuals
# estimate initial value for w
se<-matrix(0,n); # squared error
alphase<-matrix(0,n-1,k);
for (t in 2:(n+1)){ #<- time = 1 to n;
se[t-1]<-(response[t-1]-(t(F)%*%alpha_j[,t]))^2;
if (t>2){ # <- 2 to n;
past<-G%*%alpha_j[,t-1];
for (index in 1:k){
alphase[t-2,1]<-(alpha_j[(2*1)+1,t]-past[(2*1)+1])^2;
}
}
}
# variance initial value
shape<-(n/2)-1;
scale<-sum(se)/2;
vartheta = scale/(shape-1) # based on mean
# seasonal initial values
shape<-((n-1)/2)-1;
scale<-sum(alphase)/2;
w=(scale/(shape-1))/(tau[2]^2) # use mean rather than random sampling
initial.value=c(vartheta,rep(w,n.season))
return(initial.value)
}
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.