diff --git a/DESCRIPTION b/DESCRIPTION
new file mode 100644
index 0000000..7342b52
--- /dev/null
+++ b/DESCRIPTION
@@ -0,0 +1,22 @@
+Package: FracKrigingR
+Type: Package
+Title: Spatial Multivariate Data Modeling
+Version: 1.0.0
+Authors@R: c(person("Neringa", "Urbonaite", email = "neringa.urbonaite@mif.vu.lt", role = c("aut", "cre")),
+             person("Leonidas", "Sakalauskas", email = "leonidas.sakalauskas@mif.vu.lt", role = c("aut")))
+Copyright: Vilnius University Institute of Data Science and Digital
+        Technologies
+Author: Neringa Urbonaite [aut, cre],
+  Leonidas Sakalauskas [aut]
+Maintainer: Neringa Urbonaite <neringa.urbonaite@mif.vu.lt>
+Description: Aim is to provide fractional Brownian vector field generation algorithm, Hurst parameter estimation method and fractional kriging model for multivariate data modeling.
+License: GPL-2
+Encoding: UTF-8
+URL: https://github.com/NidaGreen/FracKriging
+Imports: psych, clusterGeneration, graphics, stats
+Suggests: knitr, gstat, sp, rmarkdown, raster
+RoxygenNote: 7.1.2
+NeedsCompilation: no
+Packaged: 2021-11-05 13:42:32 UTC; nerin
+Repository: CRAN
+Date/Publication: 2021-11-08 08:40:08 UTC
diff --git a/MD5 b/MD5
new file mode 100644
index 0000000..9e3c6e0
--- /dev/null
+++ b/MD5
@@ -0,0 +1,10 @@
+32294b9e9a7f24a92988e7df36ad82f5 *DESCRIPTION
+fe09b42dbe196b45447dbdda66b56dba *NAMESPACE
+b08b0433c4b7157251526e7c53186fc5 *R/FracField.R
+541d30f9ba5af5b93719d46b4ff8569e *R/FracKrig.R
+ced8e8af452d155f5223e51372e2a01d *R/FracMatrix.R
+fe7b2558ab5f66425c67e2c0037228ed *R/MaxLikelihood.R
+44fce02cf14a6696bd669d1356a8e419 *man/FracField.Rd
+ccf970a32272d9283dab2747a8056d06 *man/FracKrig.Rd
+ab6009b9ab364363bbee028ee6ca657a *man/FracMatrix.Rd
+5e5b9f3db7643cedf7d0eaefe6a36ed2 *man/MaxLikelihood.Rd
diff --git a/NAMESPACE b/NAMESPACE
new file mode 100644
index 0000000..099ba2f
--- /dev/null
+++ b/NAMESPACE
@@ -0,0 +1,8 @@
+
+
+export(FracField)
+export(FracKrig)
+export(FracMatrix)
+export(MaxLikelihood)
+importFrom("graphics", "mtext", "points")
+importFrom("stats", "rnorm")
diff --git a/R/FracField.R b/R/FracField.R
new file mode 100644
index 0000000..a5682dc
--- /dev/null
+++ b/R/FracField.R
@@ -0,0 +1,60 @@
+#' @title FracField
+#'
+#' @description Generates fractional Brownian vector field data
+#'
+#' @param  K number of observations
+#' @param  m number of criteria
+#' @param  H Hurst parameter (a real in interval [0,1))
+#' @param  X Coordinates
+
+#' @return Returns a fractional Brownian vector field matrix.
+#'
+#' @examples
+#' # Load FracKrigingR library
+#' library(FracKrigingR)
+#' # generate Coordinates
+#'    p=2; K=10;
+#'    X<-matrix(0,ncol=p, nrow=K)
+#'    for(j in 1:p){
+#'      for(i in 1:K){
+#'        X[i,j] = rnorm(1, 0, 1)
+#'      }
+#'    }
+#'    # generate fractional Brownian vector field
+#'    H = 0.5; m = 3
+#'    FracField(K,m,H,X)
+#'
+#' @importFrom graphics mtext
+#'
+#' @importFrom points stats rnorm
+#' @export
+
+
+
+
+FracField<- function(K,m,H,X){
+  suppressWarnings({
+   positive_definite_matrix <- (clusterGeneration::genPositiveDefMat(dim = m)$Sigma)
+   cholecky_matrix<- t(chol(positive_definite_matrix))
+
+   E<-as.matrix( c(rep(c(1), times = K)))
+
+    gB<-function(H){
+      gg<-((E%*%t(E))-((FracMatrix(H,K,X)*c(t(E)%*%solve(FracMatrix(H,K,X))%*%E))/2))
+      gg
+
+    }
+
+    vector<- function(m,K){
+    ksi<-matrix(0,ncol =m,nrow=K)
+    for (l in seq(from =1,to =m, by=1)){ksi[,l]<-rnorm(K,0,1)};
+    ksi
+    }
+
+    y<-t(cholecky_matrix%*%t(vector(m,K))%*%((chol(gB(H)))))
+
+    return(y)
+})
+}
+
+
diff --git a/R/FracKrig.R b/R/FracKrig.R
new file mode 100644
index 0000000..fa3d5ae
--- /dev/null
+++ b/R/FracKrig.R
@@ -0,0 +1,106 @@
+#' @title FracKrig
+#'
+#' @description Performs extrapolation for spatial multivariate data
+#'
+#' @param   X Coordinates
+#' @param   Z observations
+#' @param   Xnew Coordinates of points where the prognosis should be made
+#' @param   H Hurst parameter (a real in interval [0,1))
+#'
+#' @return   Returns a matrix of fractional kriging prognosis.
+#'
+#'
+#' @examples
+#'
+#' library(sp)
+#' library(gstat)
+#'  data(meuse)
+#'  xy<-cbind(meuse$x,meuse$y)
+
+#'  X<-xy[1:50,]
+#'  min_max_norm <- function(x) {
+#'      (x - min(x)) / (max(x) - min(x))
+#'  }
+#'  normalize <- function(x) {
+#'  return ((x - min(x)) / (max(x) - min(x)))
+#'  }
+
+
+#'  dat<-cbind(meuse[3],meuse[4],meuse[5])
+#'  data<-dat[51:100,]
+#'  zz1 <- as.data.frame(lapply(dat, normalize))
+#'  data1=as.data.frame(lapply(as.data.frame(data), normalize))
+#'  Z<-as.matrix(zz1[1:50,])
+
+# Load FracKrigingR library
+#' library(FracKrigingR)
+#'  K<-50
+
+#' #Hurst parameter estimation
+#'  H<-0.2
+#'  Xnew<-xy[51:100,]
+#'  results<- FracKrig(X,Z,Xnew,H)
+#'  denormalize <- function(x, bottom, top){
+#'     (top - bottom) * x + bottom
+#'  }
+
+#'z1 = denormalize(
+#'  results[,1], top = max(data[,1]), bottom = min(data[,1])
+#')
+#'z2 = denormalize(
+#' results[,2], top = max(data[,2]), bottom = min(data[,2])
+#')
+#'z3 = denormalize(
+#'  results[,3], top = max(data[,3]), bottom = min(data[,3])
+#')
+
+
+#'RMSE<-function(z,prognosis){
+#'  rmse<-sqrt(((1/(length(z))))*sum((z-prognosis)^2))
+#'  rmse
+#'}
+
+#'Cd<-RMSE(data[,1],z1)
+#'Cu<-RMSE(data[,2],z2)
+#'Pb<-RMSE(data[,3],z3)
+
+
+#'Cd
+#'Cu
+#'Pb
+#'
+#' @export
+#'
+
+
+FracKrig<- function(X,Z,Xnew,H){
+  X<-as.matrix(X)
+  Z<-as.matrix(Z)
+  Xnew<-as.matrix(Xnew)
+  E<-c(rep(1,K))
+  E<-as.matrix(E)
+  K<-nrow(X)
+  p<-ncol(X)
+  P<-nrow(Xnew)
+  m<-ncol(Z)
+
+  prognosis<-function(xpr,H){
+
+    b<-c()
+    for (i in seq(from =1, to=K, by=1)){
+      b<-c(b,(((xpr - t(X[i,]))%*%t(xpr-t(X[i,])))^H))
+    }
+
+    prog<-(as.matrix(t(Z)%*%solve(FracMatrix(H,K,X))))%*%((as.matrix(b))+(E%*%(((1-c(t(E)%*%solve(FracMatrix(H,K,X))%*%(as.matrix(b))))/c(t(E)%*%solve(FracMatrix(H,K,X))%*%E)))))
+    prog
+
+  }
+
+
+  final_result = matrix(c(0), nrow= P, ncol =m)
+  for (i in seq(from =1, to=P, by=1)){
+
+    final_result[i,]<-t(prognosis(Xnew[i,],H))
+  }
+  final_result
+}
diff --git a/R/FracMatrix.R b/R/FracMatrix.R
new file mode 100644
index 0000000..7f846c0
--- /dev/null
+++ b/R/FracMatrix.R
@@ -0,0 +1,37 @@
+#' @title  FracMatrix
+#'
+#' @description  Fractional distance matrix
+#'
+#' @param    H Hurst parameter (a real in interval [0,1))
+#' @param    K number of observations
+#' @param    X Coordinates
+#'
+#' @return   Returns a fractional distance matrix, which depends on the Hurst parameter.
+#'
+#' @examples
+#' # Load FracKrigingR library
+#' library(FracKrigingR)
+#' #Fractional Brownian vector field
+#'     K = 10; H = 0.5; p = 2
+#' #Generate coordinates
+#'     X<-matrix(0,ncol=p, nrow=K)
+#'     for(j in 1:p){
+#'         for(i in 1:K){
+#'             X[i,j] = rnorm(1, 0, 1)
+#'         }
+#'     }
+#'     FracMatrix(H, K, X)
+#' @export
+
+
+
+
+FracMatrix<-function(H, K, X){
+  a<-matrix(0,nrow=K,ncol=K)
+  for (i in 1:K){
+    for (j in 1:K){
+      a[i,j]<- ((t(X[i,]-X[j,]))%*%((X[i,])-X[j,]))^H
+    }
+  }
+  a
+}
diff --git a/R/MaxLikelihood.R b/R/MaxLikelihood.R
new file mode 100644
index 0000000..56718e4
--- /dev/null
+++ b/R/MaxLikelihood.R
@@ -0,0 +1,110 @@
+#' @title MaxLikelihood
+#'
+#' @description  Maximum likelihood method for Hurst parameter estimation of multivariate data
+#'
+#' @param X Coordinates
+#' @param Z Observations
+#'
+#' @return Returns the estimate of the Hurst parameter (a real in [0,1))
+#'  and a graph indicating the minimized maximum likelihood function with the Hurst parameter.
+#' @examples
+#' # Load FracKrigingR library
+#' library(FracKrigingR)
+#' # generate Coordinates
+#'    p<-2; K<-20;
+#'    X<-matrix(0,ncol=p, nrow=K)
+#'    for(j in 1:p){
+#'      for(i in 1:K){
+#'        X[i,j] = rnorm(1, 0, 1)
+#'      }
+#'    }
+#'    # generate fractional Brownian vector field
+#'    H <- 0.8; m <- 3
+#'    Z<-FracField(K,m,H,X)
+#'   # Hurst parameter estimation
+#'    MaxLikelihood(X,Z)
+#'
+#' @export
+
+MaxLikelihood<- function(X,Z){
+  Z<-as.matrix(Z)
+  K = nrow(X)
+  m = ncol(Z)
+  E<-c(rep(1,K))
+  E<-as.matrix(E)
+
+
+  golden.section.search = function(f, lower.bound, upper.bound, tolerance)
+  {
+    golden.ratio = 2/(sqrt(5) + 1)
+
+    x1 = upper.bound - golden.ratio*(upper.bound - lower.bound)
+    x2 = lower.bound + golden.ratio*(upper.bound - lower.bound)
+
+    f1 = f(x1)
+    f2 = f(x2)
+
+    ite = 0
+    while (abs(upper.bound - lower.bound) > tolerance)
+    {
+      ite = ite + 1
+      if (f2 > f1)
+      {
+        upper.bound = x2
+        x2 = x1
+        f2 = f1
+        x1 = upper.bound - golden.ratio*(upper.bound - lower.bound)
+        f1 = f(x1)
+      }
+      else
+      {
+        lower.bound = x1
+        x1 = x2
+        f1 = f2
+        x2 = lower.bound + golden.ratio*(upper.bound - lower.bound)
+        f2 = f(x2)
+      }
+    }
+    minn = (lower.bound + upper.bound)/2
+    minn
+  }
+
+
+
+
+  A<-function(H){
+
+    a<-matrix(0,nrow=K,ncol=K)
+    for (i in 1:K){
+      for (j in 1:K){
+        a[i,j]<- ((t(X[i,]-X[j,]))%*%((X[i,])-X[j,]))^H
+      }
+    }
+    a
+  }
+
+  bc<-function(H){
+    bc<-((1/K)*(((((t(Z)%*%solve(A(H))%*%E)%*%t((t(Z)%*%solve(A(H))%*%E)))/c(t(E)%*%solve(A(H))%*%E)))-(t(Z)%*%solve(A(H))%*%Z)))
+    bc
+  }
+
+  Ff<-function(H){
+    f<-(log(det(bc(H)))/m)+log(abs(det(A(H))))/K
+    f
+  }
+
+
+  t<-seq(from=0.01,to=0.9,by=0.01)
+  vect<-c()
+  for (i in t){
+    vect<-c(vect,Ff(i))
+  }
+  gsc<-golden.section.search(Ff, 0.01, 0.9, 1e-5)
+  ff<-Ff(gsc)
+
+  plot(t,vect,type="l", col="darkblue", lwd=1, xlab="", ylab="",cex.axis=1.3)
+  points(gsc,ff,pch=8,col="darkgreen",cex=2)
+  mtext(text = "Hurst Parameter", side = 1,line = 2.4,font=1,cex=1.4)
+  mtext(text = "Maximum Likelihood", side = 2, line = 2.4,font=1,cex=1.4)
+  gsc
+}
diff --git a/man/FracField.Rd b/man/FracField.Rd
new file mode 100644
index 0000000..72551c4
--- /dev/null
+++ b/man/FracField.Rd
@@ -0,0 +1,39 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/FracField.R
+\name{FracField}
+\alias{FracField}
+\title{FracField}
+\usage{
+FracField(K, m, H, X)
+}
+\arguments{
+\item{K}{number of observations}
+
+\item{m}{number of criteria}
+
+\item{H}{Hurst parameter (a real in interval [0,1))}
+
+\item{X}{Coordinates}
+}
+\value{
+Returns a fractional Brownian vector field matrix.
+}
+\description{
+Generates fractional Brownian vector field data
+}
+\examples{
+# Load FracKrigingR library
+library(FracKrigingR)
+# generate Coordinates
+   p=2; K=10;
+   X<-matrix(0,ncol=p, nrow=K)
+   for(j in 1:p){
+     for(i in 1:K){
+       X[i,j] = rnorm(1, 0, 1)
+     }
+   }
+   # generate fractional Brownian vector field
+   H = 0.5; m = 3
+   FracField(K,m,H,X)
+
+}
diff --git a/man/FracKrig.Rd b/man/FracKrig.Rd
new file mode 100644
index 0000000..8493826
--- /dev/null
+++ b/man/FracKrig.Rd
@@ -0,0 +1,71 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/FracKrig.R
+\name{FracKrig}
+\alias{FracKrig}
+\title{FracKrig}
+\usage{
+FracKrig(X, Z, Xnew, H)
+}
+\arguments{
+\item{X}{Coordinates}
+
+\item{Z}{observations}
+
+\item{Xnew}{Coordinates of points where the prognosis should be made}
+
+\item{H}{Hurst parameter (a real in interval [0,1))}
+}
+\value{
+Returns a matrix of fractional kriging prognosis.
+}
+\description{
+Performs extrapolation for spatial multivariate data
+}
+\examples{
+
+library(sp)
+library(gstat)
+ data(meuse)
+ xy<-cbind(meuse$x,meuse$y)
+ X<-xy[1:50,]
+ min_max_norm <- function(x) {
+     (x - min(x)) / (max(x) - min(x))
+ }
+ normalize <- function(x) {
+ return ((x - min(x)) / (max(x) - min(x)))
+ }
+ dat<-cbind(meuse[3],meuse[4],meuse[5])
+ data<-dat[51:100,]
+ zz1 <- as.data.frame(lapply(dat, normalize))
+ data1=as.data.frame(lapply(as.data.frame(data), normalize))
+ Z<-as.matrix(zz1[1:50,])
+library(FracKrigingR)
+ K<-50
+#Hurst parameter estimation
+ H<-0.2
+ Xnew<-xy[51:100,]
+ results<- FracKrig(X,Z,Xnew,H)
+ denormalize <- function(x, bottom, top){
+    (top - bottom) * x + bottom
+ }
+z1 = denormalize(
+ results[,1], top = max(data[,1]), bottom = min(data[,1])
+)
+z2 = denormalize(
+results[,2], top = max(data[,2]), bottom = min(data[,2])
+)
+z3 = denormalize(
+ results[,3], top = max(data[,3]), bottom = min(data[,3])
+)
+RMSE<-function(z,prognosis){
+ rmse<-sqrt(((1/(length(z))))*sum((z-prognosis)^2))
+ rmse
+}
+Cd<-RMSE(data[,1],z1)
+Cu<-RMSE(data[,2],z2)
+Pb<-RMSE(data[,3],z3)
+Cd
+Cu
+Pb
+
+}
diff --git a/man/FracMatrix.Rd b/man/FracMatrix.Rd
new file mode 100644
index 0000000..6c12cd3
--- /dev/null
+++ b/man/FracMatrix.Rd
@@ -0,0 +1,35 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/FracMatrix.R
+\name{FracMatrix}
+\alias{FracMatrix}
+\title{FracMatrix}
+\usage{
+FracMatrix(H, K, X)
+}
+\arguments{
+\item{H}{Hurst parameter (a real in interval [0,1))}
+
+\item{K}{number of observations}
+
+\item{X}{Coordinates}
+}
+\value{
+Returns a fractional distance matrix, which depends on the Hurst parameter.
+}
+\description{
+Fractional distance matrix
+}
+\examples{
+# Load FracKrigingR library
+library(FracKrigingR)
+#Fractional Brownian vector field
+    K = 10; H = 0.5; p = 2
+#Generate coordinates
+    X<-matrix(0,ncol=p, nrow=K)
+    for(j in 1:p){
+        for(i in 1:K){
+            X[i,j] = rnorm(1, 0, 1)
+        }
+    }
+    FracMatrix(H, K, X)
+}
diff --git a/man/MaxLikelihood.Rd b/man/MaxLikelihood.Rd
new file mode 100644
index 0000000..d8f3314
--- /dev/null
+++ b/man/MaxLikelihood.Rd
@@ -0,0 +1,38 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/MaxLikelihood.R
+\name{MaxLikelihood}
+\alias{MaxLikelihood}
+\title{MaxLikelihood}
+\usage{
+MaxLikelihood(X, Z)
+}
+\arguments{
+\item{X}{Coordinates}
+
+\item{Z}{Observations}
+}
+\value{
+Returns the estimate of the Hurst parameter (a real in [0,1))
+ and a graph indicating the minimized maximum likelihood function with the Hurst parameter.
+}
+\description{
+Maximum likelihood method for Hurst parameter estimation of multivariate data
+}
+\examples{
+# Load FracKrigingR library
+library(FracKrigingR)
+# generate Coordinates
+   p<-2; K<-20;
+   X<-matrix(0,ncol=p, nrow=K)
+   for(j in 1:p){
+     for(i in 1:K){
+       X[i,j] = rnorm(1, 0, 1)
+     }
+   }
+   # generate fractional Brownian vector field
+   H <- 0.8; m <- 3
+   Z<-FracField(K,m,H,X)
+  # Hurst parameter estimation
+   MaxLikelihood(X,Z)
+
+}