# Dimension of the transition matrix
dim=2
# Parameter definition
a=0.2; b=0.4;
# Definition of the transition matrix
P=matrix(c(1-a,a,b,1-b),nrow=dim,ncol=dim,byrow=TRUE)
# Number of time steps
N=100
# Encoding of chain values
Z=array(N+1);
for(ll in seq(1,N)) {
Z[1]=sample(dim,size=1,prob=P[2,])
# Random simulation of Z[j+1] given Z[j]
for (j in seq(1,N)) Z[j+1]=sample(dim,size=1,prob=P[Z[j],])
Y=array(N+1);
S=0;
# Computation of the average over the l first steps
for(l in seq(1,N+1)) { Z[l]=Z[l]-1; S=S+Z[l]; Y[l]=S/l; }
X=array(N+1); for(l in seq(1,N+1)) { X[l]=l-1; }
par(mfrow=c(2,1))
plot(X,Y,type="l",yaxt="n",xaxt="n",xlim=c(0,N),xlab="",ylim=
c(0,1),ylab="",xaxs="i",
col="black",main="",bty="n")
segments(0,a/(a+b),N,a/(a+b))
axis(2,pos=0,at=c(0.0,0.1,0.2,0.3,0.4,0.5,0.6,0.7,0.8,0.9,1.0))
axis(1,pos=0,at=seq(0,N,10),outer=TRUE)
plot(X,Z,type="o",xlab="",ylab="",xlim=c(0,N),yaxt="n",xaxt="n",xaxs="i",
col="black",main="",pch=20,bty="n")
axis(1,pos=1,at=seq(0,N+1,10),outer=TRUE,padj=-4,tcl=0.5)
axis(1,pos=0,at=seq(0,N+1,10),outer=TRUE)
axis(2,las=2,at=0:1)
readline(prompt = "Pause. Press <Enter> to continue...")
}


#MCMC
poisson.metro=function(lambda,i,n) #i:starting value
{
y=seq(n)
for(k in 1:n)
{
u1=runif(1)
j =if(u1<.5)
ifelse(i==0,i,i-1) else i+1 #symm. random walk if i>0
r =switch(i+2-j,lambda/j,1,i/lambda)
u2 =runif(1)
new=if(r>=1) j else
{if(u2<r) j else i}
i=new
y[k]=i
}
return(y)
}

yy=poisson.metro(6,100,1000)
 plot(yy,type="l")
 mean(yy)
hist(yy,prob=T,breaks=c(0:14)-0.5)
lines(c(0:14)-0.5,dpois(6,c(0:14)),col=4)

#elágazó folyamat
library(Branching)
d=1 #dimension
a=1.1 #param
n=100

Dists.i <- data.frame( name=rep( "pois", d*d ),
param1=a,
stringsAsFactors=FALSE )
rA <- rBGWM(Dists.i, "independents", d, n, z0=1) #starting pop
plot(rA$o.tt.s,type="l")


a=1
n=100

Dists.i <- data.frame( name=rep( "pois", d*d ),
param1=a,
stringsAsFactors=FALSE )
rA <- rBGWM(Dists.i, "independents", d, n, z0=1) #starting pop
plot(rA$o.tt.s,type="l")


#Brown mozgás
library(e1071)
x <- rwiener(end = 20000, frequency = 1)
 plot(x,type="l")

#iterált log tétel, nem teljesen jó
N=40000
x <- rwiener(end =N, frequency = 1)
t=c(1:N)
plot(x/t,type="l",ylim=c(-1,1),log="x")

lines(t,2*log(log(t))/sqrt(t),col=4)

lines(t,2*log(log(t))/sqrt(t),col=4)

#acf
n=1000
x=c(1:n)+rnorm(1000)
acf(f)
#nem stac.

#####
#AR(1)

n=1000
x=rep(0,times=n)
a=0.5
x[1]=rnorm(1)
for (i in 2:n)
x[i]=a*x[i-1]+(1-a^2)*rnorm(1)
plot(x,type="l")
acf(x)

n=1000
x=rep(0,times=n)
a=-0.9
x[1]=rnorm(1)
for (i in 2:n)
x[i]=a*x[i-1]+(1-a^2)*rnorm(1)
plot(x,type="l")
acf(x)

####
#AR(1) egyenletes stac. eloszlással
N=1000
x=rep(0,times=N)
for (i in 2:N)
x[i]=(x[i-1]+rbinom(1,1,0.5))*.5
plot(x[1:100],type="l")
hist(x)

####
#####
#AR(2)

n=1000
x=rep(0,times=n)
a=0.9
b=-0.2
x[1]=rnorm(1)
x[2]=rnorm(1)
for (i in 3:n)
x[i]=a*x[i-1]+b*x[i-2]+(1-a^2-b^2)*rnorm(1)
plot(x,type="l")

acf(x)
sd(x[1:500])
sd(x[501:1000])

#####
#becslések
#

#AR(1), itt minden szép és jó

a=0.9

N=100 #ism
n=1000 #eddig megyünk
ered=matrix(0,N,n)
for (i in 1:N)
{s=rep(0,times=n+1); s[1]=rnorm(1)
for (j in 1:(n))
{s[j+1]=a*s[j]+(1-a^2)*rnorm(1)
ered[i,j]=mean(s[2:(j+1)])}
}
ve=rep(0,times=n)
szor=ve
for (i in 1:n)
{ve[i]=mean(ered[,i])
szor[i]=sd(ered[,i])
}
par(mfrow=c(1,2))
plot(ve,type="l")
plot(szor,type="l")


#nem stac: véletlen bolyongás
N=100 #ism
n=100 #eddig megyünk
ered=matrix(0,N,n)
for (i in 1:N)
{s=rep(0,times=n+1); for (j in 1:n)
{s[j+1]=s[j]+round(runif(1))*2-1
ered[i,j]=mean(s[2:(j+1)])}
}
ve=rep(0,times=n)
szor=ve
for (i in 1:n)
{ve[i]=mean(ered[,i])
szor[i]=sd(ered[,i])
}
par(mfrow=c(1,2))
plot(ve,type="l")
plot(szor,type="l")
#nem konzisztens

#stac., de mégsem konzisztens (nem rövid a memória):
#trigonometrikus fv., véletlen időpontból kezdve + véletlen hiba

N=100 #ism
n=100 #eddig megyünk - 12 a periódushossz
ered=matrix(0,N,n)
for (i in 1:N)
{s=sin((c(1:12)+trunc(runif(1,0,12)))*3.14159/6)+runif(1)
s2=rep(s,times=n/12+1)
for (j in 1:n)
{
ered[i,j]=mean(s2[1:(j+1)])}
}
ve=rep(0,times=n)
szor=ve
for (ii in 1:n)
{ve[ii]=mean(ered[,ii])
szor[ii]=sd(ered[,ii])
}
par(mfrow=c(1,2))
plot(ve,type="l")
plot(szor,type="l")
#nem konzisztens

####
#tesztek a fehér zaj tulajdonságra

n=400
x=rep(0,times=n)
a=0.20
x[1]=rnorm(1)
for (i in 2:n)
x[i]=a*x[i-1]+(1-a^2)*rnorm(1)

acf(x)
lines(c(0:40),a^c(0:40),col=3)
ar(x, aic = TRUE, order.max = NULL,
   method = c("yule-walker")) #Akaike krit.

bt=Box.test(x, lag = 10, type = c("Box-Pierce"), fitdf = 0)
 bt

bt=Box.test(x, lag = 10, type = c("Ljung-Box"), fitdf = 0)

bt

#BIC a forecast csomagban van
library(forecast)

##############
#ARCH modell

library(fGarch)
spec = garchSpec(model = list(omega = 1.0e-6,
alpha = 0.10, beta = 0)) #arch(1)
x=garchSim(spec, n = 100, n.start=100)
par(mfrow=c(1,1))
plot(x,type="l")
par(mfrow=c(1,2))
acf(x)
acf(x^2)

spec = garchSpec(model = list(omega = 1.0e-6,
alpha = 0.97, beta = 0)) #arch(1)
x=garchSim(spec, n = 100, n.start=100)
par(mfrow=c(1,1))
plot(x,type="l")
par(mfrow=c(1,2))
acf(x)
acf(x^2)


spec = garchSpec(model = list(omega = 1.0e-6,
alpha = 0.5, beta = 0.48)) #garch(1,1) 
x=garchSim(spec, n = 100, n.start=100)
par(mfrow=c(1,1))
plot(x,type="l")
par(mfrow=c(1,2))
acf(x)
acf(x^2)
####################
#többdimenziós próbák
####################
library(Hotelling)
data(container.df)
fit = hotelling.test(.~gp, data = container.df)
fit

subs.df = container.df
subs.df$gp = rep(c(1,2,1,2), c(5,5,5,5))
fitPerm = hotelling.test(Al+Fe~gp, data  = subs.df, perm =  TRUE)
fitPerm
plot(fitPerm)


fitPerm = hotelling.test(Al+Fe~gp, data  = container.df, perm =  TRUE)
fitPerm
plot(fitPerm)

#######
#főkomponens analízis
#######
install.packages('factoextra',repos='http://cran.us.r-project.org')

library(factoextra) #az ábrákhoz


data(decathlon2)
decathlon2.active <- decathlon2[1:23, 1:10]
head(decathlon2.active[, 1:6]) #
res.pca <- prcomp(decathlon2.active, scale = TRUE) #átskálázzuk az eredményeket
fviz_eig(res.pca) #scree plot
eig.val <- get_eigenvalue(res.pca)
eig.val #sajátértékek

fviz_pca_var(res.pca,
             col.var = "contrib", # Color by contributions to the PC
             gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
             repel = TRUE     # Avoid text overlapping
             ) #az egyes változók (sportágak) ábrázolása
res.var <- get_pca_var(res.pca)
res.var$coord          # Coordinates
res.var$contrib        # Contributions to the PCs
res.var$cos2           # Quality of representation 

fviz_pca_ind(res.pca,
             col.ind = "cos2", # Color by the quality of representation
             gradient.cols = c("#00AFBB", "#E7B800", "#FC4E07"),
             repel = TRUE     # Avoid text overlapping
             ) #az egyes versenyzők ábrázolása

#############
#faktoranalízis
#############

v1 <- c(1,1,1,1,1,1,1,1,1,1,3,3,3,3,3,4,5,6)

v2 <- c(1,2,1,1,1,1,2,1,2,1,3,4,3,3,3,4,6,5)

v3 <- c(3,3,3,3,3,1,1,1,1,1,1,1,1,1,1,5,4,6)

v4 <- c(3,3,4,3,3,1,1,2,1,1,1,1,2,1,1,5,6,4)

v5 <- c(1,1,1,1,1,3,3,3,3,3,1,1,1,1,1,6,4,5)

v6 <- c(1,1,1,2,1,3,3,3,4,3,1,1,1,2,1,6,5,4)

m1 <- cbind(v1,v2,v3,v4,v5,v6)
 #v1~v2;v3~v4;v5~v6
cor(m1)

factanal(m1, factors = 3)


#####
#Lasso és ridge reg
#####

swiss <- datasets::swiss
x <- model.matrix(Fertility~., swiss)[,-1]
y <- swiss$Fertility
lambda <- 10^seq(10, -2, length = 100)
install.packages('glmnet',repos='http://cran.us.r-project.org')

library(glmnet)
set.seed(489)
train = sample(1:nrow(x), nrow(x)/2)
test = (-train)
ytest = y[test]
swisslm <- lm(Fertility~., data = swiss)
coef(swisslm)

ridge.mod <- glmnet(x, y, alpha = 0, lambda = lambda)
predict(ridge.mod, s = 0, type = 'coefficients')[1:6,] #szinte ugyanaz

swisslm <- lm(Fertility~., data = swiss, subset = train)
ridge.mod <- glmnet(x[train,], y[train], alpha = 0, lambda = lambda)
cv.out <- cv.glmnet(x[train,], y[train], alpha = 0) #alpha=0:ridge
bestlam <- cv.out$lambda.min
#make predictions
ridge.pred <- predict(ridge.mod, s = bestlam, newx = x[test,])
s.pred <- predict(swisslm, newdata = swiss[test,])
#check MSE
mean((s.pred-ytest)^2)
## jobb a ridge:
mean((ridge.pred-ytest)^2)
out = glmnet(x[train,],y[train],alpha = 0)
predict(ridge.mod, type = "coefficients", s = bestlam)[1:6,] #kisebb absz. értékűek

lasso.mod <- glmnet(x[train,], y[train], alpha = 1, lambda = lambda)
cv.out <- cv.glmnet(x[train,], y[train], alpha = 1) #alpha=1:lasso
bestlaml <- cv.out$lambda.min
lasso.pred <- predict(lasso.mod, s = bestlaml, newx = x[test,])
mean((lasso.pred-ytest)^2) #rosszabb, mint a ridge
lasso.coef  <- predict(lasso.mod, type = 'coefficients', s = bestlaml)[1:6,]
lasso.coef