diff --git a/src/TR.jl b/src/TR.jl
index 7d75b51..c79eb94 100644
--- a/src/TR.jl
+++ b/src/TR.jl
@@ -2,6 +2,65 @@
 
 
 
+"""
+### TO DO: ADD FRACTIONAL DERIVATIVE REGULARIZATION <-- Check that it is correctly added:) ###
+
+    regularizationMatrix(X; regType="L2", regParam1=0, regParam2=1e-14)
+    regularizationMatrix(p::Int64; regType="L2", regParam1=0, regParam2=1e-14)
+
+Calculates and returns square regularization matrix.
+
+Inputs:
+    - X/p       : Size of regularization matrix or data matrix (size of reg. mat. will then be size(X,2)
+    - regType   : "L2" (returns identity matrix)
+                  "bc" (boundary condition, forces zero on right endpoint for derivative regularization) or
+                  "legendre" (no boundary condition, but fills out reg. mat. with lower order polynomial trends to get square matrix)
+                  "std" (standardization, FILL OUT WHEN DONE)
+    - regParam1 : Int64, Indicates degree of derivative regularization (0 gives L\\_2)
+    - regParam2 : For regType=="plegendre" added polynomials are multiplied by sqrt(regParam2)
+
+Output: Square regularization matrix
+"""
+function regularizationMatrix(X; regType="L2", regParam1=0, regParam2=1e-14)
+
+if regType == "std"
+    regParam2 = Diagonal(vec(std(X, dims=1)));
+end
+
+regMat = regularizationMatrix(size(X,2); regType, regParam1, regParam2)
+
+return regMat
+end
+
+function regularizationMatrix(p::Int64; regType="L2", regParam1=0, regParam2=1e-14)
+    
+if regType == "bc" # Discrete derivative with boundary conditions
+    regMat  = [I(p); zeros(regParam1,p)]; for i = 1:regParam1 regMat = diff(regMat, dims = 1); end
+elseif regType == "legendre" # Fill in polynomials in bottom row(s) to get square matrix
+    regMat  = [I(p); zeros(regParam1,p)]; for i = 1:regParam1 regMat = diff(regMat, dims = 1); end
+    P, _    = plegendre(regParam1-1, p);   
+    regMat[end-regParam1+1:end,:] = sqrt(regParam2) * P;
+elseif regType == "L2"
+    regMat = I(p);
+elseif regType == "std" # Standardization
+    regMat = regParam2;
+elseif regType == "GL" # Grünwald-Letnikov fractional derivative regulariztion
+    # regParam1 is alpha (order of fractional derivative)
+    C = ones(p)*1.0;
+    for k in 2:p
+        C[k] = (1-(regParam1+1)/(k-1)) * C[k-1];
+    end
+
+    regMat = zeros(p,p);
+
+    for i in 1:p
+        regMat[i:end, i] = C[1:end-i+1];
+    end
+end
+
+return regMat
+end
+
 """
     function TRLooCV