diff --git a/DESCRIPTION b/DESCRIPTION
new file mode 100644
index 0000000..020d898
--- /dev/null
+++ b/DESCRIPTION
@@ -0,0 +1,32 @@
+Package: flap
+Type: Package
+Title: Forecast Linear Augmented Projection
+Version: 0.1.0
+Authors@R: 
+      person(given = "Yangzhuoran Fin", 
+             family = "Yang", 
+             email = "yangyangzhuoran@gmail.com", 
+             role = c("aut", "cre"), 
+             comment = c(ORCID = "0000-0002-1232-8017")) 
+Description: The Forecast Linear Augmented Projection (flap) method reduces 
+    forecast variance by adjusting the forecasts of multivariate time series to 
+    be consistent with the forecasts of linear combinations (components) of the 
+    series by projecting all forecasts onto the space where the linear 
+    constraints are satisfied. The forecast variance can be reduced 
+    monotonically by including more components. For a given number of 
+    components, the flap method achieves maximum forecast variance reduction 
+    among linear projections. 
+License: GPL (>= 3)
+Encoding: UTF-8
+URL: https://github.com/FinYang/flap
+BugReports: https://github.com/FinYang/flap/issues/
+Imports: corpcor, utils
+RoxygenNote: 7.2.3
+Suggests: forecast, stats
+NeedsCompilation: no
+Packaged: 2024-02-27 22:36:38 UTC; fyan0012
+Author: Yangzhuoran Fin Yang [aut, cre]
+    (<https://orcid.org/0000-0002-1232-8017>)
+Maintainer: Yangzhuoran Fin Yang <yangyangzhuoran@gmail.com>
+Repository: CRAN
+Date/Publication: 2024-02-28 19:10:09 UTC
diff --git a/MD5 b/MD5
new file mode 100644
index 0000000..6bdabe5
--- /dev/null
+++ b/MD5
@@ -0,0 +1,8 @@
+4421ac36e00b36d917ca599bc5fe2cd7 *DESCRIPTION
+dac34b3e67f836b8b4bc57be733fba1e *NAMESPACE
+e1c6577add934633543df00339084a49 *NEWS.md
+4d69bd05792836eb9a458ff1baad4b3c *R/flap.R
+ca46a024684c12e7bfe4dddd4b2ee422 *R/projection.R
+dca4700102c39c5bae07e6c467f5b4c4 *R/utils.R
+283114137f3cf4e0099825ce5cf49034 *README.md
+73b2e08cd8b4376d9d858ec1f206226d *man/flap.Rd
diff --git a/NAMESPACE b/NAMESPACE
new file mode 100644
index 0000000..765d9f8
--- /dev/null
+++ b/NAMESPACE
@@ -0,0 +1,5 @@
+# Generated by roxygen2: do not edit by hand
+
+S3method(as.data.frame,flap)
+S3method(print,flap)
+export(flap)
diff --git a/NEWS.md b/NEWS.md
new file mode 100644
index 0000000..4c8945c
--- /dev/null
+++ b/NEWS.md
@@ -0,0 +1,3 @@
+# flap 0.1.0
+
+* Initial CRAN submission.
diff --git a/R/flap.R b/R/flap.R
new file mode 100644
index 0000000..64f417e
--- /dev/null
+++ b/R/flap.R
@@ -0,0 +1,107 @@
+#' Forecast Linear Augmented Projection
+#'
+#' Reduces forecast variance by adjusting the forecasts of multivariate time
+#' series to be consistent with the forecasts of linear combinations (components)
+#' of the series by projecting all forecasts onto the space where the linear
+#' constraints are satisfied.
+#'
+#' @param fc An \eqn{h} by \eqn{m} matrix of base forecasts of the original
+#' series to be projected. \eqn{h} is the forecast horizon and \eqn{m} is the
+#' total number of series.
+#' @param fc_comp An \eqn{h} by \eqn{p} matrix of base forecasts of the components
+#' used in the projection. \eqn{h} is the forecast horizon and \eqn{p} is the
+#' total number of components.
+#' @param Phi A \eqn{p} by \eqn{m} weight matrix mapping the original series into
+#' the components such that \eqn{c_t = \Phi z_t} where \eqn{c_t} is the vector of
+#' components and \eqn{z_t} is the vector of original series.
+#' @param res A \eqn{T} by \eqn{m} (in-sample) forecast residual matrix of the
+#' original series.
+#' @param res_comp A \eqn{T} by \eqn{p} (in-sample) forecast residual matrix of
+#' the components.
+#' @param p The number of components to use in the projection. The default is
+#' trying all the possible number of components capped at the number provided in
+#' the forecast.
+#'
+#' @return A list of class \code{flap} with each element containing a \eqn{h} by
+#' \eqn{m} matrix of projected forecast of the original series for the corresponding
+#' number of components \code{p}.
+#'
+#' @examples
+#' # Generate example data
+#' # T = 70, m = 20
+#' train <- matrix(rnorm(70 * 20),ncol = 20)
+#'
+#' # Obtain the forecast and the residual of the original series
+#' mdl <- apply(train, 2, forecast::ets)
+#' fc <- vapply(mdl, function(mdl) forecast::forecast(mdl, h=12)$mean,
+#'              FUN.VALUE = numeric(12))
+#' res <- vapply(mdl, residuals, FUN.VALUE = numeric(70))
+#'
+#' # Obtain components and their forecasts and residuals
+#' pca <- stats::prcomp(train, center = FALSE, scale. = FALSE)
+#' mdl_comp <- apply(pca$x, 2, forecast::ets)
+#' fc_comp <- vapply(mdl_comp, function(mdl) forecast::forecast(mdl, h=12)$mean,
+#'                   FUN.VALUE = numeric(12))
+#' res_comp <- vapply(mdl_comp, residuals,
+#'                    FUN.VALUE = numeric(nrow(pca$x)))
+#' Phi <- t(pca$rotation)
+#'
+#' # flap!
+#' flap(fc, fc_comp, Phi, res, res_comp)
+#'
+#' @export
+flap <- function(fc, fc_comp, Phi, res, res_comp,
+                 p = seq_len(ncol(fc_comp))) {
+  W <- get_W(res, res_comp, p)
+  proj_fc <- project(
+    cbind(fc,
+          fc_comp),
+    W = W,
+    Phi = Phi,
+    p = p)
+  class(proj_fc) <- c("flap", class(proj_fc))
+  proj_fc
+}
+
+#' @export
+as.data.frame.flap <- function(x, row.names = NULL, optional = FALSE, ...){
+  mat <- do.call(rbind, x)
+  df <- as.data.frame(mat)
+  df$p <- rep(as.integer(names(x)), each = nrow(x[[1]]))
+  df$h <- rep(seq_len(nrow(x[[1]])), times = length(x))
+  df
+}
+
+#' @export
+print.flap <- function(x, ...) {
+  cat("Forecast Linear Augmented Projection\n")
+  cat("A named list of numeric matrices of projected forecasts\n")
+  cat("------------")
+
+  cs <- "Num. of Series:"
+  ns <- nchar(cs)
+  cc <- "Num. of Components:"
+  nc <- nchar(cc)
+  ch <- "Num. of Forecast Horizons:"
+  nh <- nchar(ch)
+
+  cat("\n", paste0(cs, strrep(" ", nh-ns)), "m =", ncol(x[[1]]))
+  cat("\n", paste0(cc, strrep(" ", nh-nc)), "p = ")
+  print_consecutive(as.integer(names(x)))
+  cat("\n", paste0(ch), nrow(x[[1]]))
+  cat("\n------------\n")
+
+  utils::str(x, vec.len = 2, give.attr = FALSE, list.len = 5)
+  invisible(x)
+}
+
+print_consecutive <- function(nums) {
+  g <- cumsum(c(1L, diff(nums) != 1))
+  r <- rle(g)
+  end <- cumsum(r$lengths)
+  start <- c(1L, 1L + end[-length(end)])
+  out <- paste(nums[start], nums[end],sep = "-")
+  out[start == end] <- nums[start[start == end]]
+  cat(out, sep = ", ")
+  invisible(nums)
+}
diff --git a/R/projection.R b/R/projection.R
new file mode 100644
index 0000000..77e981c
--- /dev/null
+++ b/R/projection.R
@@ -0,0 +1,37 @@
+project <- function(fc, W, Phi, p) {
+  C_all <- cbind(-Phi, diag(nrow(Phi)))
+  m <- ncol(fc) - nrow(Phi)
+  proj_fc <- lapply(
+    asplit(fc, 1), function(fc){
+      mapply(function(p, W){
+        C <- block(C_all, p, m+p)
+        WtC <- tcrossprod(W, C)
+        bf <- c(fc[seq_len(m+p)])
+        (bf -tcrossprod(WtC, t(solve(C %*% WtC, C))) %*% bf)[seq_len(m),]
+      },
+      p = p,
+      W = W,
+      SIMPLIFY = FALSE)
+    })
+
+  proj_fc <- lapply(proj_fc, function(x) do.call(cbind, x))
+  proj_fc <- list2array(proj_fc)
+  proj_fc <- aperm(proj_fc, c(3, 1, 2))
+  colnames(proj_fc) <- colnames(fc)[seq_len(m)]
+  proj_fc <- array2list(proj_fc)
+  names(proj_fc) <- p
+  proj_fc
+}
+
+block <- function(mat, m, n = m){
+  mat[seq_len(m), seq_len(n), drop = FALSE]
+}
+
+get_W <- function(res_ori, res_com, p) {
+  m <- NCOL(res_ori)
+  res <- cbind(res_ori, res_com)
+  res <- res[!apply(res, 1, anyNA),]
+  lapply(
+    p,
+    function(pp) corpcor::cov.shrink(res[,seq_len(m+pp)], verbose = FALSE))
+}
diff --git a/R/utils.R b/R/utils.R
new file mode 100644
index 0000000..0df144a
--- /dev/null
+++ b/R/utils.R
@@ -0,0 +1,15 @@
+list2array <- function(xlist){
+  d1 <- unique(vapply(xlist, NROW, numeric(1)))
+  if(length(d1) != 1) stop("Different row number")
+  d2 <- unique(vapply(xlist, NCOL, numeric(1)))
+  if(length(d2) != 1) stop("Different col number")
+  array(unlist(xlist), dim = (c(d1, d2, length(xlist))))
+}
+
+array2list <- function(xarray){
+  out <- vector("list", length= dim(xarray)[[3]])
+  for(i in seq_len(dim(xarray)[[3]])) {
+    out[[i]] <- xarray[,,i]
+  }
+  out
+}
diff --git a/README.md b/README.md
new file mode 100644
index 0000000..b1d6961
--- /dev/null
+++ b/README.md
@@ -0,0 +1,99 @@
+
+<!-- README.md is generated from README.Rmd. Please edit that file -->
+
+# flap
+
+<!-- badges: start -->
+
+[![CRAN
+status](https://www.r-pkg.org/badges/version/flap)](https://CRAN.R-project.org/package=flap)
+[![R-CMD-check](https://github.com/FinYang/flap/actions/workflows/R-CMD-check.yaml/badge.svg)](https://github.com/FinYang/flap/actions/workflows/R-CMD-check.yaml)
+[![Licence](https://img.shields.io/badge/licence-GPL--3-blue.svg)](https://www.gnu.org/licenses/gpl-3.0.en.html)
+<!-- badges: end -->
+
+The goal of `flap` is to provide the Forecast Linear Augmented
+Projection method that can reduce forecast variance by adjusting the
+forecasts of multivariate time series to be consistent with the
+forecasts of linear combinations (components) of the series by
+projecting all forecasts onto the space where the linear constraints are
+satisfied.
+
+## Installation
+
+You can install the **stable** version from
+[CRAN](https://cran.r-project.org/package=flap).
+
+``` r
+install.packages("flap")
+```
+
+You can install the **development** version from
+[Github](https://github.com/FinYang/flap)
+
+``` r
+# install.packages("remotes")
+remotes::install_github("FinYang/flap")
+```
+
+## Example
+
+This is a basic workflow to flap:
+
+``` r
+## The following pacakges are required to run this example
+# install.packages("tidyr")
+# install.packages("ggplot2")
+# install.packages("forecast")
+# install.packages("fpp2")
+
+library(flap)
+library(tidyr)
+library(ggplot2)
+
+# Obtain the forecast and the residual of the original series
+mdl <- apply(fpp2::visnights, 2, forecast::ets)
+#> Registered S3 method overwritten by 'quantmod':
+#>   method            from
+#>   as.zoo.data.frame zoo
+fc <- vapply(mdl, function(mdl) forecast::forecast(mdl, h=12)$mean, 
+             FUN.VALUE = numeric(12))
+res <- vapply(mdl, residuals, 
+              FUN.VALUE = numeric(nrow(fpp2::visnights)))
+
+# Obtain components and their forecasts and residuals
+pca <- stats::prcomp(fpp2::visnights, center = FALSE, scale. = FALSE)
+mdl_comp <- apply(pca$x, 2, forecast::ets)
+fc_comp <- vapply(mdl_comp, function(mdl) forecast::forecast(mdl, h=12)$mean, 
+                  FUN.VALUE = numeric(12))
+res_comp <- vapply(mdl_comp, residuals, 
+                   FUN.VALUE = numeric(nrow(pca$x)))
+Phi <- t(pca$rotation)
+
+# flap!
+proj_fc <- flap(fc, fc_comp, Phi, res, res_comp)
+proj_fc
+#> Forecast Linear Augmented Projection
+#> A named list of numeric matrices of projected forecasts
+#> ------------
+#>  Num. of Series:            m = 20
+#>  Num. of Components:        p = 1-20
+#>  Num. of Forecast Horizons: 12
+#> ------------
+#> List of 20
+#>  $ 1 : num [1:12, 1:20] 7.8 7.91 ...
+#>  $ 2 : num [1:12, 1:20] 7.64 7.76 ...
+#>  $ 3 : num [1:12, 1:20] 7.64 7.78 ...
+#>  $ 4 : num [1:12, 1:20] 7.39 7.48 ...
+#>  $ 5 : num [1:12, 1:20] 7.39 7.49 ...
+#>   [list output truncated]
+
+# Plot
+if(interactive()) {
+  proj_fc %>% 
+    as.data.frame() %>% 
+    pivot_longer(!c(h, p)) %>% 
+    ggplot(aes(x = h, y = value, colour = p, group = p)) +
+    geom_line() +
+    facet_wrap("name", scales = "free_y")
+}
+```
diff --git a/man/flap.Rd b/man/flap.Rd
new file mode 100644
index 0000000..96f48c3
--- /dev/null
+++ b/man/flap.Rd
@@ -0,0 +1,66 @@
+% Generated by roxygen2: do not edit by hand
+% Please edit documentation in R/flap.R
+\name{flap}
+\alias{flap}
+\title{Forecast Linear Augmented Projection}
+\usage{
+flap(fc, fc_comp, Phi, res, res_comp, p = seq_len(ncol(fc_comp)))
+}
+\arguments{
+\item{fc}{An \eqn{h} by \eqn{m} matrix of base forecasts of the original
+series to be projected. \eqn{h} is the forecast horizon and \eqn{m} is the
+total number of series.}
+
+\item{fc_comp}{An \eqn{h} by \eqn{p} matrix of base forecasts of the components
+used in the projection. \eqn{h} is the forecast horizon and \eqn{p} is the
+total number of components.}
+
+\item{Phi}{A \eqn{p} by \eqn{m} weight matrix mapping the original series into
+the components such that \eqn{c_t = \Phi z_t} where \eqn{c_t} is the vector of
+components and \eqn{z_t} is the vector of original series.}
+
+\item{res}{A \eqn{T} by \eqn{m} (in-sample) forecast residual matrix of the
+original series.}
+
+\item{res_comp}{A \eqn{T} by \eqn{p} (in-sample) forecast residual matrix of
+the components.}
+
+\item{p}{The number of components to use in the projection. The default is
+trying all the possible number of components capped at the number provided in
+the forecast.}
+}
+\value{
+A list of class \code{flap} with each element containing a \eqn{h} by
+\eqn{m} matrix of projected forecast of the original series for the corresponding
+number of components \code{p}.
+}
+\description{
+Reduces forecast variance by adjusting the forecasts of multivariate time
+series to be consistent with the forecasts of linear combinations (components)
+of the series by projecting all forecasts onto the space where the linear
+constraints are satisfied.
+}
+\examples{
+# Generate example data
+# T = 70, m = 20
+train <- matrix(rnorm(70 * 20),ncol = 20)
+
+# Obtain the forecast and the residual of the original series
+mdl <- apply(train, 2, forecast::ets)
+fc <- vapply(mdl, function(mdl) forecast::forecast(mdl, h=12)$mean,
+             FUN.VALUE = numeric(12))
+res <- vapply(mdl, residuals, FUN.VALUE = numeric(70))
+
+# Obtain components and their forecasts and residuals
+pca <- stats::prcomp(train, center = FALSE, scale. = FALSE)
+mdl_comp <- apply(pca$x, 2, forecast::ets)
+fc_comp <- vapply(mdl_comp, function(mdl) forecast::forecast(mdl, h=12)$mean,
+                  FUN.VALUE = numeric(12))
+res_comp <- vapply(mdl_comp, residuals,
+                   FUN.VALUE = numeric(nrow(pca$x)))
+Phi <- t(pca$rotation)
+
+# flap!
+flap(fc, fc_comp, Phi, res, res_comp)
+
+}