Giovanni Ferretti 1.0 2026-06-23 Functions to work with tensors C10 pdfdoc:TenTS.pdf examples examples/trade_data.gdt examples/fama_french.gdt examples/application_to_fama_french.inp /* -doc~ This function transform a matrix into a tensor of dimensions specified by the matrix "dims". It invokes tsz_special() if the tensor has 3 dimensions or less, otherwise it starts a recursion. The order of observations is given by the vec() of original matrix INPUT: -table: original matrix -dims: dimensions of resulting tensor OUTPUT: -ret: bundle containing the tensor */ ### general case: meaning that the tensor order is greater than 3 n = nelem(dims) bundle ret = empty d = vec(dims)' if prodr(dims)!= nelem(X) funcerr "Number of observations is not equal to product of Dimensions!" endif # if tensor order is not greater than 3, call the tsz_special function if n < 4 ret = tsz_special(X, dims) else particular_case = n>=4 && d[n-1]==1 ? 1 : 0 # rename size of n-th dimension and n-1 dimension k = dims[n-1] h = dims[n] ini = 1 fin = k # create an empty array with size equal to the size of n-th dimension arrays ret.tensor = array(h) # consider sizes up to n-1 dimension dims2 = dims[1:n-1] # Consider a temporary tensor with n-1 dimensions instead of n: last dimension will be the product # between n-1 and n-th one dims2[n-1] *= dims[n] tmp = tensorize(X, dims2) if particular_case tmp = very_particular( tmp, dims2) endif # start a loop on the n-th dimension and populate the tensor that will be returned # with the values from the temporary tensor (that has n-1 dimensions) organized in # n dimensions if nelem(tmp.dims) != nelem(dims) loop i = 1 .. h ret.tensor[i] = tmp.tensor[ini:fin] ini = fin + 1 fin += k endloop endif endif # add dims and time_id as metadata ret.dims = dims if exists(time_id) ret.time_id = time_id endif # add strings for dimension names if exists(dimnames) ret.dimnames = dimnames endif return ret /* -doc~ The following function extracts the dimensions of a tensor INPUT: -a: bundle containing the tensor OUTPUT: -ret: matrix containing dimensions of the tensor */ bundle tmp = a tmp.dim = {} tmp_elem = tmp.tensor type = typename(tmp_elem) tmp.type = type if type != "matrix" loop i = 1.. 100 tmp = go_deep(tmp) type = tmp.type if type == "matrix" break endif endloop endif matrix ret = tmp.dim tmp_elem2 = tmp.tensor ret ~= {cols(tmp_elem2), rows(tmp_elem2)} ret = mreverse(ret,1) return ret /* -doc~ Transform a tensor of N dimensions to a tensor of 3 dimensions, where the third dimension is given by the product of all dimensions greater or equal than 3 of the input tensor INPUT: tensor_bundle: bundle containing the tensor OUTPUT: list_mat: object of type matrices where third dimension size equal to the product of all dimensions greater or equal than 3 */ matrix dims = exists(tensor_bundle.dims)?tensor_bundle.dims:get_tensor_dims(tensor_bundle) scalar N = cols(dims) matrices list_mat = array(0) arr = tensor_bundle.tensor # first take care of dimensions greater than 5 if N>=5 loop i=N..5 --decr # define empty array of arrays arrays temp = array(0) loop n=1..nelem(arr) # pupulate the array of arrays, flattening dimensions n arrays temp = temp+arr[n] endloop # update the array of array with the flattened one arrays arr = temp endloop # Now populate the list of matrices loop i=1..nelem(arr) list_mat = list_mat + arr[i] endloop elif N==4 # if the dimension of the tensor is 4, skip the first step loop i=1..nelem(arr) z = arr[i] list_mat = list_mat + arr[i] endloop # if the dimension of the tensor is 3, directly return the list of matrices elif N==3 list_mat = arr else funcerr "Tensor has less than 3 dimensions!" endif return list_mat /* -doc~ INPUT: tensor_bundle: bundle containing the tensor to be vectorized OUTPUT: vec_ten: the vectorized tensor */ # extract tensor dimensions matrix dims = exists(tensor_bundle.dims)?tensor_bundle.dims:get_tensor_dims(tensor_bundle) scalar N = cols(dims) # transfom the tensor into a 3d tensor, with third dim equal to product of dimensons # greater than 3 if N>=3 matrices list_mat = tensor_3d(tensor_bundle) ret = vec(flatten(list_mat,0)) else ret = vec(tensor_bundle.tensor) endif return ret /* -doc~ The following function unfolds an N-th order tensor into its n-th mode. It is more efficient with respect to kolda_bader_unfolding. note that if "tensor" is a matrices object, this is equivalent to flatten(tensor, mode) INPUT tensor_bundle: bundle containing the tensor n: scalar specifying the dimension to unfold OUTPUT ret: matrix wit unfolded tensor */ matrix dims = exists(tensor_bundle.dims)?tensor_bundle.dims:get_tensor_dims(tensor_bundle) scalar N = cols(dims) if N > 2 matrices tensor3 = tensor_3d(tensor_bundle) if sum(dims.=1)>0 && n>2 ret = kolda_bader_unfolding( compute_idmat(dims), dims, vec(flatten(tensor3,0)) ,n) elif n==1 ret = flatten(tensor3) elif n==2 ret = flatten(tensor3,1)' elif n==3 # create empty matrix to store results matrix ret = {} # consider a number of iterations # equal to the product of dimensions from 4 to N scalar iterations = nelem(tensor3)/dims[n] scalar start_count = 1 scalar end_count = 0 # At each iteration, consider a number of matrices equal to # the third dimension, flatten them and horizontally stack them loop j = 1 .. iterations end_count += dims[n] ret ~= flatten(tensor3[start_count:end_count],2)' start_count = end_count+1 endloop else # create empty matrix to store results matrix ret = {} # define useful quantities scalar start_count = 1 scalar end_count = 0 # At each iteration, consider a number of matrices # equal to the product of dimensions from 3 to n-1 scalar roll_count = prodr(dims[3:n-1]) # For dimension n, number of iterations is the product of dimensions # from n+1 to N scalar iterations = n<N?prodr(dims[n+1:N]):1 loop j=1 .. iterations matrix temp = {} loop i=1 .. dims[n] end_count += roll_count temp |= vec(flatten(tensor3[start_count:end_count]))' start_count = end_count+1 endloop ret ~= temp endloop endif else # handle case where N <= 2 ret = n==1?tensor_bundle.tensor:tensor_bundle.tensor' endif return ret /* -doc~ This function changes the dimension order of a tensor INPUT: tensor_bundle: bundle containing the tensor new_order: matrix containing the new order of dimensions OUTPUT: tensor_bundle: tensor with updated dimensions order */ # get actual dimensions of the tensor matrix old_dims = exists(tensor_bundle.dims)?tensor_bundle.dims:get_tensor_dims(tensor_bundle) errorif(cols(new_order)!=cols(old_dims),"Number of dimensions is not the same") #compute index matrix for each obserrvation matrix id_mat = compute_idmat(old_dims) #vectorize the tensor matrix vec_ten = vectorize_tensor(tensor_bundle) # change the column-order for the index matrix matrix new_id_mat = id_mat[,new_order] matrix new_dims = old_dims[,new_order] #compute first dimension mode with the new order matrix mode_1 = kolda_bader_unfolding(new_id_mat, new_dims, vec_ten, 1) #tensorize the mode bundle ret = tensorize(mode_1,new_dims) ret.dims = new_dims # fix order for dimension names if exists(tensor_bundle.dimnames) ret += _(dimnames = permute_dimnames(tensor_bundle.dimnames, new_order) ) endif if exists(tensor_bundle.time_id) check = dsort( ({tensor_bundle.time_id}.=new_order).*seq(1,nelem(old_dims)) ) ret.time_id = check[1] endif return ret /* -doc~ This function implements the product between a tensor And an array of matrices. It uses property described in "Tensor Decomposition and Applications" (Kolda and bader 2009, pag. 461) INPUT tensor_bundle: bundle containing the tensor loads: list of loading matrices transpose: if the loadings need to be transposed OUTPUT reduced_tensor: tensor of new dimensions */ # Extract tensor dimensions matrix dims = exists(tensor_bundle.dims)?tensor_bundle.dims:get_tensor_dims(tensor_bundle) scalar N = cols(dims) # check if all dimensions have a loading matrix, if not, add an Identity matrix of the respective size loop i = 1 .. nelem(loads) if nelem(loads[i])==0 loads[i] = I(dims[i]) endif endloop # Unfold the tensor along first dimension matrix Xn = unfold(tensor_bundle,1) # Extract dimensions of the tensor after multiplication # If transpose is true, transpose the loading matrices if transpose loads = transp_array(loads) endif matrix reduced_dims = zeros(1,N) loop i=1..N reduced_dims[,i] = rows(loads[i]) endloop # Do the product using aformentioned property matrix load_prod = 1 loop i= nelem(loads) .. 2 --decr load_prod = load_prod**loads[i] endloop matrix Yn = loads[1]*Xn*load_prod' # Split back into tensor after multiplication ret = tensorize(Yn,reduced_dims) if exists(tensor_bundle.time_id) ret.time_id = tensor_bundle.time_id endif return ret /* -doc~ Computes averages of a tensor of order N along dimension specified in index INPUT tensor_bundle: bundle containing the tensor index: scalar denoting dimension to consider for average computation squeeze: 1 if the dimension used for average computation is to be dropped OUTPUT ret: tensor of order N (N-1 if squeeze=1) containing the averages along dimension specified in index */ matrix dims = exists(tensor_bundle.dims)?tensor_bundle.dims:get_tensor_dims(tensor_bundle) scalar N = cols(dims) if N==3 X = tensor_bundle.tensor r = dims[1] c = dims[2] n = dims[3] if index == 1 averages = mshape(meanc(flatten(X)), c, n) elif index == 2 averages = mshape(meanr(flatten(X, 1)), r, n) elif index == 3 averages = mshape(meanr(flatten(X, 2)), r, c) endif else matrix averages = meanc(unfold(tensor_bundle,index)) #bundle ret = tensorize(avg_mode, dims[-index]) endif # return either an N-1 order tensor, or a tensor with dimension used for mean eq. to 1 if squeeze dims = dims[-index] else dims[index] = 1 endif bundle ret = tensorize(averages,dims) if exists(tensor_bundle.dimnames) && squeeze dimnames = tensor_bundle.dimnames[-index] ret += _(dimnames = dimnames) elif exists(tensor_bundle.dimnames) dimnames = tensor_bundle.dimnames dimnames[index] = defarray("") ret += _(dimnames = dimnames) endif # add time_id if present if exists(tensor_bundle.time_id) && (squeeze==0 || tensor_bundle.time_id != index) ret.time_id = tensor_bundle.time_id endif return ret /* -doc~ Computes sum of a tensor of order N along dimension specified in index INPUT tensor_bundle: bundle containing the tensor index: scalar denoting dimension to consider for sum computation squeeze: 1 if the dimension used for sum is to be dropped OUTPUT ret: tensor of order N (N-1 if squeeze=1) containing the sum along dimension specified in index */ matrix dims = exists(tensor_bundle.dims)?tensor_bundle.dims:get_tensor_dims(tensor_bundle) scalar N = cols(dims) if N==3 X = tensor_bundle.tensor r = dims[1] c = dims[2] n = dims[3] if index == 1 averages = mshape(sumc(flatten(X)), c, n) elif index == 2 averages = mshape(sumr(flatten(X, 1)), r, n) elif index == 3 averages = mshape(sumr(flatten(X, 2)), r, c) endif else matrix averages = sumc(unfold(tensor_bundle,index)) #bundle ret = tensorize(avg_mode, dims[-index]) endif # return either an N-1 order tensor, or a tensor with dimension used for mean eq. to 1 if squeeze dims = dims[-index] else dims[index] = 1 endif bundle ret = tensorize(averages,dims) if exists(tensor_bundle.dimnames) && squeeze dimnames = tensor_bundle.dimnames[-index] ret += _(dimnames = dimnames) elif exists(tensor_bundle.dimnames) dimnames = tensor_bundle.dimnames dimnames[index] = defarray("") ret += _(dimnames = dimnames) endif # add time_id if present if exists(tensor_bundle.time_id) && (squeeze==0 || tensor_bundle.time_id != index) ret.time_id = tensor_bundle.time_id endif return ret /* -doc~ Computes standard deviation of a tensor of order N along dimension specified in index INPUT tensor_bundle: bundle containing the tensor index: scalar denoting dimension to consider for standard deviation computation squeeze: 1 if the dimension used for standard deviation is to be dropped OUTPUT ret: tensor of order N (N-1 if squeeze=1) containing the standard deviation along dimension specified in index */ matrix dims = exists(tensor_bundle.dims)?tensor_bundle.dims:get_tensor_dims(tensor_bundle) scalar N = cols(dims) if N==3 X = tensor_bundle.tensor r = dims[1] c = dims[2] n = dims[3] if index == 1 sd_mode = mshape(sdc(flatten(X)), c, n) elif index == 2 sd_mode = mshape(sdc(flatten(X, 1)'), r, n) elif index == 3 sd_mode = mshape(sdc(flatten(X, 2)'), r, c) endif else matrix sd_mode = sdc(unfold(tensor_bundle,index)) endif # return either an N-1 order tensor, or a tensor with dimension used for mean eq. to 1 if squeeze dims = dims[-index] else dims[index] = 1 endif # tensorize and add dimnames if they are there bundle ret = tensorize(sd_mode,dims) if exists(tensor_bundle.dimnames) && squeeze dimnames = tensor_bundle.dimnames[-index] ret += _(dimnames = dimnames) elif exists(tensor_bundle.dimnames) dimnames = tensor_bundle.dimnames dimnames[index] = defarray("") ret += _(dimnames = dimnames) endif # add time_id if present if exists(tensor_bundle.time_id) && (squeeze==0 || tensor_bundle.time_id != index) ret.time_id = tensor_bundle.time_id endif return ret /* -doc~ Print all matrices of the tensor labelling them with respective dimensions from 3 to N INPUT tensor_bundle: bundle containing the tensor OUTPUT void */ matrix dims = exists(tensor_bundle.dims)?tensor_bundle.dims:get_tensor_dims(tensor_bundle) matrix id_mat = compute_idmat(dims) scalar N = cols(dims) if N >= 3 # transform the N dimensional tensor into a 3d tensor matrices list_mat = tensor_3d(tensor_bundle) # print matrices string name = argname(tensor_bundle) scalar counter = 0 loop i=1..nelem(list_mat) counter += dims[1]*dims[2] matrix dims_i = id_mat[counter,] matrix dims_i = dims_i[3:N] string label = build_label(dims_i) printf "%s%s\n\n", name, label printf "%10.4g", list_mat[i] printf "\n\n" endloop else print tensor_bundle.tensor endif /* -doc~ Computes elementwise sum between two tensors of order N INPUT A: bundle containing the first tensor B: bundle containing the second tensor OUTPUT ret: tensor of order N containing the elementwise sum between A and B */ matrix dims = exists(A.dims)?A.dims:get_tensor_dims(A) matrix a = vectorize_tensor(A) matrix b = vectorize_tensor(B) matrix sum_ab = a+b bundle ret = tensorize(sum_ab, dims) if exists(A.time_id) ret.time_id = A.time_id endif if exists(A.dimnames) && !which_name ret += _(dimnames = A.dimnames) elif exists(B.dimnames) && which_name ret += _(dimnames = B.dimnames) endif return ret /* -doc~ Computes elementwise difference between two tensors of order N INPUT A: bundle containing the first tensor B: bundle containing the second tensor OUTPUT ret: tensor of order N containing the elementwise difference between A and B */ matrix dims = exists(A.dims)?A.dims:get_tensor_dims(A) matrix a = vectorize_tensor(A) matrix b = vectorize_tensor(B) matrix diff_ab = a-b bundle ret = tensorize(diff_ab, dims) if exists(A.time_id) ret.time_id = A.time_id endif if exists(A.dimnames) && !which_name ret += _(dimnames = A.dimnames) elif exists(B.dimnames) && which_name ret += _(dimnames = B.dimnames) endif return ret /* -doc~ Computes elementwise product between two tensors of order N INPUT A: bundle containing the first tensor B: bundle containing the second tensor OUTPUT ret: tensor of order N containing the elementwise product between A and B */ matrix dims = exists(A.dims)?A.dims:get_tensor_dims(A) matrix a = vectorize_tensor(A) matrix b = vectorize_tensor(B) matrix c = a.*b bundle ret = tensorize(c,dims) if exists(A.time_id) ret.time_id = A.time_id endif if exists(A.dimnames) && !which_name ret += _(dimnames = A.dimnames) elif exists(B.dimnames) && which_name ret += _(dimnames = B.dimnames) endif return ret /* -doc~ Computes elementwise ratio between two tensors of order N INPUT A: bundle containing the first tensor B: bundle containing the second tensor OUTPUT ret: tensor of order N containing the elementwise ratio between A and B */ matrix dims = exists(A.dims)?A.dims:get_tensor_dims(A) matrix a = vectorize_tensor(A) matrix b = vectorize_tensor(B) matrix c = a./b bundle ret = tensorize(c,dims) if exists(A.time_id) ret.time_id = A.time_id endif if exists(A.dimnames) && !which_name ret += _(dimnames = A.dimnames) elif exists(B.dimnames) && which_name ret += _(dimnames = B.dimnames) endif return ret /* -doc~ Computes outerproduct between two tensors, respectively of order N1 and N2 INPUT A: bundle containing the first tensor B: bundle containing the second tensor OUTPUT ret: tensor of order N1+N2 containing the outerproduct between A and B */ matrix dims_a = exists(A.dims)?A.dims:get_tensor_dims(A) matrix dims_b = exists(B.dims)?B.dims:get_tensor_dims(B) matrix ret_dims = dims_a~dims_b matrix a = vectorize_tensor(A) matrix b = vectorize_tensor(B) matrix c = b**a ret = tensorize(c,ret_dims) return ret /* -doc~ This function elementwise applies function specified in "op" to the two tensors A and B INPUT A: bundle containing the first tensor B: bundle containing the second tensor OUTPUT ret: tensor of order N containing the output of "op" */ matrix dims = exists(A.dims)?A.dims:get_tensor_dims(A) matrix a = vectorize_tensor(A) matrix b = vectorize_tensor(B) ret = tensorize(a @op b, dims) if exists(A.time_id) ret.time_id = A.time_id endif if exists(A.dimnames) && !which_name ret += _(dimnames = A.dimnames) elif exists(B.dimnames) && which_name ret += _(dimnames = B.dimnames) endif return ret /* -doc~ INPUT A: bundle containing the tensor to filter filter_dims: bundle containing for each dimension the ids to keep OUTPUT tensor with filtered dimensions */ matrix dims = exists(A.dims) ? A.dims : get_tensor_dims(A) scalar N = cols(dims) # check if filters make sense strings filts = getkeys(filter_dims) loop i = 1 .. nelem(filts) string n = filts[i] if @n>N funcerr "A dimension greater than the tensor order was specified!" endif endloop matrix a = vectorize_tensor(A) matrix id_mat = compute_idmat(dims) # consider only dimensions to filter (those where filter matrices are not null) matrix to_filter = {} loop i = 1 ..N if inbundle(filter_dims, "$i") to_filter |= i endif endloop # We have a total of k dimensions to filter scalar k = rows(to_filter) # Loop across the dimensions to filter: if the observation matches all the dimensions to filter, matrix check will store its row # number in the vec of original tensor check = seq(1,prodr(dims))' loop i = 1 .. k scalar dim = to_filter[i] string dim_str = sprintf("%d",dim) check = check.*sumr(id_mat[,dim] .= vec(filter_dims[dim_str])') endloop check = sort(check) zero_count = sum(check.=0) check = check[zero_count+1:rows(check),] # Compute dimensions of the filtered tensor: # if dimension was not filtered, use original dimension (from dims) # if dimension was filtered, use the filtered dimension matrix new_dims = {} loop i = 1 ..N elems = inbundle(filter_dims, "$i") if elems == 0 new_dims ~= dims[i] else new_dims ~= nelem(filter_dims["$i"]) endif endloop # if resulting tensor is a scalar, fix the dimensions scalar prod_dims = prodr(new_dims) if prod_dims==1 new_dims = {1} endif # filter out observations that need to be discarded and tensorize the others a = a[check] ret = tensorize(a, new_dims) # if dimnames are present, filter out names for discarded indices in each dimension if exists(A.dimnames) arrays new_names = filter_names(A.dimnames,new_dims, filter_dims ) ret += _(dimnames = new_names) endif # if time_id is present, keep it if exists(A.time_id) ret.time_id = A.time_id endif return ret /* -doc~ INPUT: tensors: Array of bundles containing the tensors coeffs: matrix containing coefficients for the linear combination OUTPUT: tensor resulting from the linear combination of the arrays of bundles passed as input */ bundle A = tensors[1] matrix dims = exists(A.dims)?A.dims:get_tensor_dims(A) # compute linear combination of all elements of the array of bundles matrix ret_vec = zeros(prodr(dims),1) loop i = 1 .. nelem(tensors) ret_vec += coeffs[i].*vectorize_tensor(tensors[i]) endloop # tensorize the linear combination and add meta-data bundle ret = tensorize(ret_vec, dims) if exists(A.time_id) ret.time_id = A.time_id endif if exists(A.dimnames) ret += _(dimnames = A.dimnames) endif return ret /* -doc~ Function to compute lags of a tensor along dimension specified by the bundle-key "time_id" INPUT: A: bundle containing the tensor to compute lags P: scalar denoting the order of the lag trim_zeros: scalar denoting how many time periods need to be discarded at the beginning of the sample OUTPUT: lagged tensor */ matrix dims = exists(A.dims)?A.dims:get_tensor_dims(A) scalar N = cols(dims) if exists(A.time_id) == 0 funcerr "Dimension index for time was not specified!" endif scalar time_id = A.time_id if exists(A.dimnames) names = A.dimnames endif if N > 2 # now set time as first dimension of the tensor ord = seq(1,N) new_ord = {time_id}~ord[-time_id] new_dims = dims[new_ord] A = permute(A,new_ord) matrix time_mat = mlag( unfold(A,1), P, m) # if trim_zeros is not zero remove first P obs. in the time dimensions if trim_zeros time_mat = time_mat[P+1:,] new_dims[1] = rows(time_mat) endif # tensorize and permute back to original order ret = permute_back( tensorize(time_mat, new_dims, 1), new_ord) elif N<=2 && time_id ==1 ret = trim_zeros ? _(tensor = mlag(A.tensor, P)[P+1:,] ) : _(tensor = mlag(A.tensor, P, m)) ret.time_id = time_id elif N<=2 && time_id ==2 ret = trim_zeros ? _(tensor = mlag(A.tensor', P)[P+1:,] ) : _(tensor = mlag(A.tensor', P,m) ) ret.time_id = time_id endif ret.dims = get_tensor_dims(ret) if exists(A.dimnames) && trim_zeros names[ret.time_id] = names[ret.time_id][-seq(1,P)] ret.dimnames = names elif exists(A.dimnames) ret.dimnames = A.dimnames endif return ret /* -doc~ Function to compute first differences of a tensor along dimension specified by the bundle-key "time_id" INPUT: A: bundle containing the tensor to compute first differences OUTPUT: tensors with first differences */ matrix dims = exists(A.dims)?A.dims : get_tensor_dims(A) scalar N = cols(dims) bundle ret = tensor_subtraction(A, tensor_lag(A,lag,0, NA) ) if trim_zeros ret = tensor_drill(ret, defbundle(sprintf("%d",A.time_id) , seq(lag+1, dims[A.time_id]) )) endif ret.time_id = A.time_id if exists(A.dimnames) ret += _(dimnames = A.dimnames) endif return ret /* -doc~ INPUT: A: bundle containing the tensor time series dim: dimension for the tensor unfolding filter_matrices: bundle containing dimensions to filter as keys, inedexes to retain as values aggr_mat: Dimensions to compute averages or sum on, (depending on do_means) time_mat: If the plotted matrix is to be stored do_means: 1 if averages of dimensions specified in aggr_mat are to be computed, 0 if sum OUTPUT: void: the function plots along the dim dimension */ matrix dims = exists(A.dims)?A.dims:get_tensor_dims(A) # Ensure all dimension have a value for their index if exists(A.dimnames) A.dimnames = check_missing_names(A.dimnames, dims) endif # if filter matrices, filter only the selected ids for each dimension of the tensor if exists(filter_matrices) A = tensor_drill(A,filter_matrices) endif # if do_means, compute the mean across selected dimensions, otherwise the sum if exists(aggr_mat) loop i = cols(aggr_mat) .. 1 --decr A = do_means?tensor_mean(A,aggr_mat[i]):tensor_aggr(A,aggr_mat[i]) endloop endif # unfold the tensor on the time dimension: resulting matrix has time on rows matrix time_mat = unfold(A,dim) if exists(A.dimnames) cnameset(time_mat, plot_names(A, dim) ) xtics = A.dimnames[dim] endif # plot data along selected dimension if exists(xtics) plot time_mat options time-series with-lines fit=none single-yaxis literal set key outside right literal set xtics rotate by 45 right printf "set xtics (%s)", xtics_lab(xtics) #literal unset ytics end plot --output=display else gnuplot --matrix=time_mat --with-lines --time-series --single-yaxis --output=display endif /* -doc~ This Function has exactly the same syntax of bwrite() */ return bwrite(A, "@path",export) /* -doc~ This Function has exactly the same syntax of bread() */ return bread("@path", import) /* -doc~ The function estimates a tensor autoregression of order P, following method described in Xi and Liao,2021 "Multi-linear Tensor Autoregressive Models". Supported estimation methods are Least squares and Maximum-likelihood. Tensor can be centered by using the boolean "demean". */ if !exists(A.time_id) funcerr "No dimension for time was specified!" endif scalar time_id = A.time_id matrix dims = exists(A.dims) ? A.dims : get_tensor_dims(A) scalar K = cols(dims) if A.time_id != K funcerr "Time dimension must be the last one!" endif # Center the tensor if demean=1 if demean demeanss = center_data(A) A = demeanss[1] means = demeanss[2] means.tensor = means.tensor[1:dims[K]-P] endif printf "Selected number of lags: %d\n", P if !MLE printf "Estimation Method: LSE\n" bundle mod = TenAr_LSE(&A, P, verbose, tol, max_iter) else printf "Estimation Method: MLE\n" bundle mod = TenAr_MLE(&A, P, verbose, tol, max_iter) endif # re-add the mean to estimates if it was removed if demean bundle ret = post_estimation(&A, mod, means) else bundle ret = post_estimation(&A, mod) endif # print results eval print_results(ret) return ret /* -doc~ This function transforms a matrix into a tensor of dimensions specified by the matrix "dims" with a maximum number of dimensions equal to 3. The order of observations is given by the vec() of original matrix INPUT: table: original matrix dims: dimensions of resulting tensor OUTPUT: tensor_bundle: bundle containing the tensor */ ### special cases: # if dimension number is 1 returns a bundle containing a vector under the key "tensor" n = nelem(dims) bundle ret = empty if n == 1 matrix ret.tensor = mshape(X, dims[1], 1) elif n == 2 # if dimension number is 2 returns a bundle containing a matrix undser the key "tensor" matrix ret.tensor = mshape(X, dims[1], dims[2]) elif n == 3 # if dimension number is 3 returns a bundle containing an array of matrices under the key "tensor" scalar r = dims[1] scalar c = dims[2] matrices ret.tensor = msplitby(mshape(X, r, c * dims[3]), c, 1) else funcerr "not a special case!" endif return ret /* -doc~ This function is invoked only in the case where in tensorize() n>=4 and the n-1 dimension has size of 1. For example, if n=4, and third dimension has size 1, doing ret.tensor[i] = tmp.tensor[ini:fin] in tensorize returns an error, since we are assigning a matrix to an element of type arrays of arrays. Following function solves this issue by using an element of type matrices as an intermediary */ scalar n = cols(dims2) scalar dn = dims2[n] arrays mat_to_array = array(dn) if n == 3 loop i = 1 .. dn matrices mat_to_array[i] = array(1) mat_to_array[i][1] = tmp.tensor[i] endloop else loop i = 1 .. dn arrays mat_to_array[i] = array(1) mat_to_array[i][1] = tmp.tensor[i] endloop endif return _(tensor = mat_to_array, dims=dims2 ) bundle ret = empty ret.dim = a.dim~{nelem(a.tensor)} if a.type =="matrices" tmp_elem = a.tensor matrix mat = typeof(tmp_elem)== 3 ? tmp_elem : tmp_elem[1] ret.tensor = mat ret.type = "matrix" elif a.type == "arrays" tmp_elem = a.tensor[1] ret.tensor = tmp_elem ret.type = typename(tmp_elem) endif return ret /* -doc~ INPUT: dims: row vector containing dimensions of the tensor OUTPUT: id_mat: matrix with number of rows equal to the total number of observations, number of columns equal to the order of the tensor. Each i-th row contains the dimensions of the i-th observation */ scalar N = cols(dims) scalar total = prodr(dims) # Call the selection id matrix "id_mat" matrix id_mat = zeros(total, N) # Loop over each mode to build the correct Kronecker-like index pattern loop i = 1..N scalar before = 1 # repeat each number this many times scalar after = 1 # number of blocks loop j=1..i-1 before*= dims[j] endloop loop j=i+1..N after*=dims[j] endloop matrix base = seq(1, dims[i])' matrix one_col = ones(before,1) matrix repeated = base ** one_col # Repeat each number 'before' times matrix tiled = ones(after,1) ** repeated # Tile the result 'after' times id_mat[,i] = tiled endloop return id_mat /* -doc~ Computes n-th mode of the tensor using algorithm from Kolda and Bader (2009), pag. 460 INPUT: id_mat: matrix with number of rows equal to the total number of observations, number of columns equal to the order of the tensor. Each i-th row contains the dimensions of the i-th observation dims: row vector containing dimensions of the tensor vec_ten: the vectorized tensor mode: dimension to unfold OUTPUT: nmode: n-th mode matrix */ scalar N = cols(dims) scalar tot_obs = rows(id_mat) r = dims[mode] c = prodr(dims[-mode]) matrix nmode = zeros(r, c) loop i=1..tot_obs matrix id = id_mat[i,] # --- use algorithm from Kolda and Bader (2009), pag. 460 scalar col = 0 loop k = 1..N if k == mode # skip the unfolding mode continue endif # J_k = product of sizes of modes < k, excluding n scalar Jk = 1 loop u=1..k-1 if u!= mode Jk*= dims[u] endif endloop col += (id[k] - 1) * Jk endloop scalar col += 1 nmode[id[mode], col] = vec_ten[i] # store i-th observation in the matrix endloop return nmode /* -doc~ This function is used in permute() in order to fix order for dimension names in the permuted tensor INPUT: dimmnames: Array of arrays containing dimension names OUTPUT: new_names: Permuted Array of arrays */ arrays new_names = array(nelem(dimnames)) loop i = 1 .. nelem(dimnames) j = new_order[i] new_names[i] = dimnames[j] endloop return new_names /* -doc~ If a permutation was done, This function permutes the tensor to the original dimension order INPUT: A: bundle containing the tensor new_order: matrix containing the permutation that was performed OUTPUT: ret: tensor with updated dimensions order */ matrix dims = exists(A.dims)?A.dims:get_tensor_dims(A) N = cols(dims) matrix invs = zeros(1,N) loop i = 1..N scalar k = new_ord[i] invs[k] = i endloop ret = permute(A,invs) return ret /* -doc~ INPUT: X: Array of matrices OUTPUT: X: Array of matrices, containing the transposed of the input one */ # In principle, one could use msplitby(flatten(X)', v), but # the efficiency gain is dubious scalar N = nelem(X) loop i=1..N X[i] = X[i]' endloop return X /* -doc~ function used in tensor_print() for printing dimensions of each matrix of the tensor */ string label = "" loop j=1..cols(dims_i) label += sprintf("[%d]", dims_i[j]) endloop return label /* -doc~ The folllowing function is used in tensor_drill() in order to filter out names of dimension indices that are filtered out by tensor_drill() INPUT dimnames: Arrays of arrays containing dimension names OUTPUT new_names: Arrays of arrays containing names of the filtered dimensions */ # order of the original tensor scalar N = cols(dims) arrays new_names = array(N) # select names that are preserved in the tensor_drill() function loop i = 1 .. N if inbundle(filter_dims, "$i") == 0 strings new_names[i] = nelem(dimnames[i])>1?dimnames[i][1:nelem(dimnames[i])]: defarray("") else strings new_names[i] = nelem(filter_dims["$i"])>1? dimnames[i][filter_dims["$i"]] : defarray(dimnames[i][filter_dims["$i"]]) endif endloop return new_names /* -doc~ This function is used in plot_names() to ensure that there are not missing values in names to plot. In practice, we substitute them with blank spaces INPUT dimnames: array of arrays containing dimension names dims: matrix containing dimensions of a tensor OUTPUT dimnames: array of arrays containing dimension names with blank spaces instead of missing names */ scalar N = cols(dims) loop i = 1 .. N if nelem(dimnames[i]) != dims[i] strings set_names = array(0) loop j = 1 .. dims[i] set_names += exists(dimnames[i][j])?dimnames[i][j] : "" endloop dimnames[i] = set_names[1:nelem(set_names)] endif endloop return dimnames /* -doc~ This function is used in tensor_plot() in order to create array of strings containing names of plotted series, that will be displayed in the legend. INPUT: A: bundle containing the tensor dim: index of the dimension that will be on the x-axis OUTPUT: s: array of strings containing names of dimensions displayed in the legend */ matrix dims = exists(A.dims)?A.dims:get_tensor_dims(A) scalar N = cols(dims) names = A.dimnames # start from first dimension that has size greater than one, and that is not the dimension that # will be on the x-axis matrix check = seq(1,N) check = sort( check.*(check.!=dim) ) tot_zero = sum(check.=0) check = check[tot_zero+1:cols(check)] strings plot_names = nelem(names[check[1]])>1? names[check[1]][1:dims[check[1]]] : names[check[1]][1] strings tmp = array(0) # then concatenate the other names check = check[-1] loop i = 1 .. nelem(check) n = check[i] strings names_dim_i = names[n][1:dims[n]] loop name = 1 .. dims[n] tmp += plot_names ~ " " ~ names_dim_i[name] endloop plot_names = tmp strings tmp = array(0) endloop return plot_names /* -doc~ This function is used in tensor_plot() to create xtics labels INPUT: lab: array of strings containing names of the dimension in the x-axis OUTPUT: s: string containing xtics and their location */ if nelem(lab)>8 scalar non_empty_freq = round(nelem(lab)/8) else non_empty_freq = 1 endif string s = "" loop i = 1 .. nelem(lab) if i%non_empty_freq ==0 s~= " '"~lab[i]~"' $i," endif endloop return s[-nelem(s)] /* -doc~ This function centers the tensor and retrieves two bundles. The first one is the centered tensor, the second one is a tensor with same dimensions of the original one, but with repeated mean per each time period. (useful if one wants to re-add the means with the function tensor_sum() ). */ matrix dims = exists(A.dims) ? A.dims : get_tensor_dims(A) scalar N = cols(dims) scalar time_size = dims[N] bundle avgs = tensor_mean(A,N, 1) matrix vecmeans = ones(time_size,1)**vectorize_tensor(avgs) bundle bigmeans = tensorize(vecmeans,dims) bundle norm_data = tensor_subtraction(A,bigmeans) bundles ret ret += norm_data # Retrieve a tensor of the same order of A with repeated means for each value of time ret += bigmeans return ret /* -doc~ This function initializes paramater matrices for a tensor autoregression of order P. Output is an array of arrays with P elements (one per lag), each containing an array of N diagonal matrices, where N is the number of dimensions of the tensor excluding time. Each matrix is square with dimensions ( dims[i] x dims[i] ), i = 1 .. N. */ scalar N = cols(dims) # block to initialize loading arrays ret = array(P) loop p = 1 .. P matrices ret[p] = array(N-1) loop i = 1 .. N-1 ret[p][i] = 0.5*I(dims[i],dims[i]) endloop endloop return ret /* -doc~ This function unfolds the original tensor in all its dimensions. Output is an array of arrays, with one element per each time period. Each element contains an array of N matrices where N is the number of tensor dimensions excluding time. Each matrix contains the matricization of the tensor at time t, along i-th dimension, t = 1 .. T, i = 1 .. N. */ matrix dims = exists(A.dims) ? A.dims : get_tensor_dims(A) scalar N = cols(dims) scalar time_id = A.time_id scalar T = dims[time_id] arrays ret = array(T) bundle xt = _(dims = dims[-N]) loop t = 1 .. T xt.tensor = A.tensor[t] matrices ret[t] = array(N-1) loop i = 1 .. N-1 ret[t][i] = unfold(xt,i) endloop endloop return ret /* -doc~ This function subtracts values predicted by other lags. */ scalar P = nelem(loads) scalar T = nelem(unfolded_tensor) scalar K = nelem(unfolded_tensor[1]) # initialize matrix of zeros of the seletced size matrix temp = zeros(rows(unfolded_tensor[sel_lag][sel_dim]),cols(unfolded_tensor[sel_lag][sel_dim]) ) loop iter_lag = 1 .. P # Consider the unfolded tensor iter_lag periods before matrix p_mode = unfolded_tensor[t-iter_lag][sel_dim] # sum all part of the tensor predicted by lags different from "sel_lag" if iter_lag!= sel_lag temp_load = loads[iter_lag][-sel_dim] # if n. of dimensions is gretar than 2, do kronecker prod. if K>2 matrix temp_load_prod = 1 loop j = nelem(temp_load) .. 1 --decr temp_load_prod = temp_load_prod**temp_load[j] endloop else temp_load_prod = temp_load endif w = p_mode*temp_load_prod' temp += loads[iter_lag][sel_dim]*w endif endloop # subtract part predicted by other lags from the observed unfolded tensor return unfolded_tensor[t][sel_dim]-temp /* -doc~ This Function updates loadings with the Least-Squares method. It modifies the values in loads, that is an array of arrays with P elements. Each of this P elements contains an array of matrices of dimension N, where N is the number of the tensor dimensions excluding time */ scalar K = cols(dims) scalar T = dims[K] scalar P = nelem(loads) loop p = 1 .. P loop i = K-1 .. 1 --decr # initialize useful quantities matrix xw = zeros(dims[i],dims[i]) matrix ww = zeros(dims[i],dims[i]) # consider all loadings of r-th rank except for the i-th dimension matrices lambda_no_d = loads[p][-i] matrix lambda_prod = {1} loop j = nelem(lambda_no_d) .. 1 --decr matrix lambda_prod = lambda_prod**lambda_no_d[j] endloop loop t = P+1 .. T mode_lag = unfolded_ten[t-p][i] if P >1 mode = subtract_mode(t, unfolded_ten, loads,p,i) else mode = unfolded_ten[t][i] endif w = mode_lag*lambda_prod' ww += w*w' xw += mode*w' endloop # update matrix of loads for dimension i temp_load = mols(xw',ww')' # divide by frob. norm for identification if it is not # last dimension if i < K-1 loads[p][i] = temp_load./sqrt(tr(temp_load'temp_load)) else loads[p][i] = temp_load endif endloop endloop return loads /* -doc~ This function computes predicted values by summing the tensor to matrix products of p-lagged tensors and loadings referring to lag p, p = 1 .. P. */ scalar P = nelem(loads) scalar N = nelem(loads[1]) matrix dims = exists(A.dims) ? A.dims : get_tensor_dims(A) scalar T = dims[N+1] # add identity matrix as loading for time dimension # (this is to avoid looping through time) loop p = 1 .. P loads[p] += I(T-P) endloop dims[N+1] = T-P bundle Ahat = tensorize(zeros(prodr(dims),1),dims) loop p = 1 .. P # select the p-lagged tensor Alag = _(tensor = A.tensor[P-p+1:T-p]) # select loadings for the p-lag matrices loads_p = loads[p][1:N+1] # sum values predicted by all P lags Ahat = tensor_sum( Ahat,tensor_times_matrices(Alag,loads_p) ) endloop return Ahat /* -doc~ This Function checks for convergence: Return max of frobenius norm of the difference between old and updated loadings */ scalar P = nelem(loads) scalar N = nelem(loads[1]) matrix norm_gains = {} loop p = 1 .. P loop i = 1 .. N matrix diff_load = loads[p][i] - old_loads[p][i] norm_gains |= sqrt(tr( diff_load'diff_load )) endloop endloop return maxc(norm_gains) /* -doc~ This Function estimates a tensor autoregression of order P with Least-Squares Estimator. Returns a bundle containing: 1) a tensor with estimated values; 2) an array of arrays with parameter matrices; 3) The unfolded original tensor (useful in post-estimation() function ). */ if !exists(A.time_id) funcerr "No dimension for time was specified!" endif scalar time_id = A.time_id matrix dims = exists(A.dims) ? A.dims : get_tensor_dims(A) scalar N = cols(dims) arrays unfolded_ten = compute_modes(&A) # intialize the loadings loads = init_loads(dims, P) # initialize the alg scalar iter = 0 scalar gain = 1 scalar converged = 0 loop while iter<max_iter && !converged iter++ old_loads = loads loads = update_loads(loads, unfolded_ten, A.dims) scalar gain = check_convergence(loads, old_loads) scalar converged = gain<= tol if verbose flush printf "Iter n. %d, gain: %g\n", iter, gain flush endif endloop # print results if converged printf "Converged at iteration %g\n", iter else printf "Maximum number of %g iterations reached. Convergence not achieved.\n", max_iter endif # store predicted values bundle Ahat = sum_pred_lags(A, loads) bundle mod = _(estimates = Ahat, param = loads, unfolded_tensor = unfolded_ten) return mod /* -doc~ This function initializes covariance matrices for a tensor autoregression of order P estimated with ML method. Output is an array of N identity matrices, where N is the number of dimensions of the tensor excluding time. */ scalar N = cols(dims) # block to initialize loading matrices ret = array(N-1) loop i = 1 .. N-1 ret[i] = I(dims[i],dims[i]) endloop return ret /* -doc~ This function updates parameter matrix for lag p, dimension i through ML method. p = 1, .. ,P and i = 1, .., N. It also updates the matrix inv_covs by computing the kroncecker product of the inverse of the array of matrices containing covariances in all the N dimensions */ scalar K = cols(dims) scalar T = dims[K] scalar P = nelem(loads) unfolded_ten = unfolded_tens[1][1:T] unfolded_res = unfolded_tens[2][1:T-P] # consider all loadings of r-th rank except for the i-th dimension matrices lambda_no_d = loads[p][-i] matrix lambda_prod = {1} # do the same for covs matrices covs_no_d = covs[-i] matrix inv_covs = {1} # compute kronecker product of selected loadings and covs loop j = nelem(lambda_no_d) .. 1 --decr lambda_prod = lambda_prod**lambda_no_d[j] inv_covs = inv_covs**ginv(covs_no_d[j],10e^-9) endloop # initialize useful quantities matrix xw = zeros(dims[i],dims[i]) matrix ww = zeros(dims[i],dims[i]) loop t = P+1 .. T mode_lag = unfolded_ten[t-p][i] if P >1 mode = subtract_mode(t, unfolded_ten, loads,p,i) else mode = unfolded_ten[t][i] endif # part for updating loadings xw += mode * inv_covs *lambda_prod*mode_lag' # xw += mode_lag*lambda_prod'*inv_covs*lambda_prod*mode_lag' ww += qform(mode_lag*lambda_prod',inv_covs) endloop # update matrix of loads for dimension i matrix ret = mols(xw',ww')' return ret /* -doc~ This function updates covariance matrix for dimension i through ML method. */ scalar K = cols(dims) scalar T = dims[K] matrix no_t = dims[-K] # product of all dimensions except for time scalar prod_no_t = prodr(no_t[-i] ) # n. of lags considered scalar P = nelem(loads) unfolded_res = unfolded_tens[2][1:T-P] # initialize useful quantities matrix s = zeros(dims[i],dims[i]) loop t = P + 1 .. T # take the i-th mode of the lagged residuals r_mode = unfolded_res[t-P][i] # Part for Updating covariances (see section 3.2 in the paper) # s+= r_mode*inv_covs*r_mode' s += qform(r_mode, inv_covs) endloop # update matrix of covs for dimension i covs[i] = s./( (T-P) *prod_no_t) /* -doc~ This function returns an array of bundles containing p-lagged tensors, p = 1, .. , P */ matrix dims = exists(A.dims) ? A.dims : get_tensor_dims(A) scalar T = dims[A.time_id] bundles ret = array(P) loop p = 1 .. P # select the p-lagged tensor ret[p] = _(tensor = A.tensor[P-p+1:T-p]) endloop return ret /* -doc~ This function sum values predicted by each of p-lagged tensors multiplied for the parameter matrices that refer to the corresponding lag. Returns the predicted tensor of the tensor autoregression of order P. */ scalar P = nelem(predicted_lags) Ahat = predicted_lags[1] loop p = 2 .. P Ahat = tensor_sum(Ahat,predicted_lags[p]) endloop return Ahat /* -doc~ This Function normalizes loadings and covariances for the tenar estimated through MLE to match results from the r package TensorTS */ scalar K = nelem(loads[1]) matrix load_norms = zeros(K-1, 1) matrix cov_norms = zeros(K-1,1) loop i = 1 .. K-1 # normalize loadings tmp_load = loads[p][i] sing_vals = svd(tmp_load) load_norms[i] = sing_vals[1,1] loads[p][i] = tmp_load./load_norms[i] # normalize covariances cov_norms[i] = mreverse(eigensym(covs[i]))[1] covs[i] = covs[i]./cov_norms[i] endloop loads[p][K] = loads[p][K].*prodc(load_norms) covs[K] = covs[K].*prodc(cov_norms) return loads arrays loads = old_loads scalar P = nelem(loads) scalar K = nelem(loads[1]) loop p = 1 .. P matrix norms = zeros(K-1,1) loop k = 1 .. K-1 tmp_load = loads[p][k] norms[k] = sqrt(tr(tmp_load'tmp_load)) loads[p][k] = loads[p][k]./norms[k] endloop loads[p][K] *= prodc(norms) endloop return loads /* -doc~ This Function estimates a tensor autoregression of order P with Maximum-Likelihood Estimator. Returns a bundle containing: 1) an array of arrays with parameter matrices; 2) a tensor with estimated values; 3) The unfolded original tensor (useful in post-estimation() function ); 4) An array of matrices containing covariances */ matrix dims = exists(A.dims) ? A.dims : get_tensor_dims(A) scalar time_id = A.time_id scalar T = dims[time_id] scalar K = cols(dims) # intitialize loadings and covariances arrays loads = init_loads(dims, P) matrices covs = init_covs(dims) # bundles containing lagged values of the tensor bundles lagged_tens = store_lagged_tens(A, P) bundles predicted_lags = array(P) arrays temp_loads = loads loop p = 1 .. P # add time loads to avoid loop through time temp_loads[p] += I(T-P) # store values predicted by each lag predicted_lags[p] = tensor_times_matrices(lagged_tens[p], temp_loads[p][1:K]) endloop # original tensor discarding first P observations bundle orig = _(tensor = A.tensor[ P+1 : T ], time_id = time_id) # final prediction is the sum of predictions from each lag bundle estimates = sum_predictions(predicted_lags) # compute residuals bundle R = tensor_subtraction(orig, estimates) # unfold original tensor and residual in all periods and dimensions bundles to_unfold = defarray(A, R) arrays unfolded_tens = array(2) loop i = 1 .. 2 unfolded_tens[i] = compute_modes(&to_unfold[i]) endloop # initialize values scalar converged = 0 scalar gain = 1 scalar iter = 0 matrix inv_covs = {} flush # Estimation loop while !converged && iter < max_iter iter++ old_loads = loads loop p = 1 .. P loop i = K-1 .. 1 --decr # UPDATE LOADINGS loads[p][i] = update_loads_ml( loads, covs, unfolded_tens, dims,i, p, &inv_covs) temp_loads[p][i] = loads[p][i] # update value predicted by p-th lag predicted_lags[p] = tensor_times_matrices(lagged_tens[p], temp_loads[p][1:K]) # sum predictions by each of the P lags bundle Ahat = P>1 ? sum_predictions(predicted_lags) : predicted_lags[1] # compute residuals bundle R = tensor_subtraction(orig, Ahat) # unfold residuals unfolded_tens[2] = compute_modes(&R) # UPDATE COVARIANCES update_covs_ml(loads, &covs, unfolded_tens, dims, i, inv_covs) endloop loads = normalize_loads_covs(loads, p, &covs) endloop # check for convergence scalar gain = check_convergence(loads, old_loads) scalar converged = gain <= delta if verbose flush printf "Iteration n. %d, gain is: %g \n", iter, gain flush endif endloop # ensure that all coefficient matrices have frob. norm of 1 except for the last dimension loads = final_norm(loads) # Post-estimation part printf "\n\n----------------------------------------------\n\n" if converged printf "Converged at iteration %g\n", iter else printf "Maximum number of %g iterations reached. Convergence not achieved.\n", max_iter endif printf "\n\n----------------------------------------------\n\n" flush bundle mod = _(param = loads, estimates = Ahat, unfolded_tensor = unfolded_tens[1], covs = covs) return mod matrix dims = E.dims scalar T = dims[cols(dims)] scalar K = cols(dims)-1 matrix dims_no_t = dims[1:K] # Initialize covariances as I(k), k = 1 .. K matrices covs = init_covs(dims) # Unfold residuals in all dimensions matrices unfolded_e = array(K) loop i = 1 .. K unfolded_e[i] = unfold(E,i) endloop # Initialize inverse of covariances matrices inv_covs = covs # flip-flop loop i = 1 .. 10 loop k = 1 .. K inv_sig = I(T,T) loop j = K .. 1 --decr if j != k inv_sig = inv_sig**inv_covs[j] endif endloop e_white = unfolded_e[k]*inv_sig Q = (e_white*unfolded_e[k]')/prodr(dims) # update covs and inv_covs covs[k] = dims[k]*Q inv_covs[k] = inv(covs[k]) endloop endloop return covs matrix dims = E.dims scalar T = dims[cols(dims)] scalar K = cols(dims)-1 dims_no_t = dims[1:K] matrices inv_covs = array(K) scalar det_penalty = 0 loop k = 1 .. K inv_covs[k] = inv(covs[k]) det_penalty += ( prodr(dims_no_t)/dims[k] ) * log(det(covs[k])) endloop # add an Identity for time to avoid loop and perform multilinear product inv_covs += I(T) whiten_E = tensor_times_matrices(E, inv_covs) # compute likelihood e_white = vectorize_tensor(whiten_E) e = vectorize_tensor(E) llik = e_white'e + T*det_penalty return -llik # P number of lags scalar P = nelem(loads) # K number of dimensions scalar K = nelem(loads[1]) matrix comp = {} # create upper part of the companion matrix loop p = 1 .. P matrix tmp = {1} loop k = K .. 1 --decr tmp = tmp**loads[p][k] endloop comp ~= tmp endloop #create bottom part of the companion matrix if P>=2 if P == 2 scalar D = prodr(dims) matrix iden = I(D) bot_comp = iden ~ zeros(D,D) comp |= bot_comp elif P>2 scalar D = prodr(dims) matrix iden = I(D) bot_comp = iden loop p = 2 .. P-1 bot_comp = diagcat(bot_comp,iden) endloop bot_comp ~= zeros((P-1)*D,D ) comp |= bot_comp endif return max(abs(eigen(comp))) /* -doc~ Given output from either TenAr_LSE or TenAr_MLE, this Function computes covariance matrix of residuals, standard error of the estimator, bic and loglikelihood. */ matrix dims = exists(A.dims) ? A.dims : get_tensor_dims(A) scalar K = cols(dims) scalar T = dims[K] loads = mod.param scalar P = nelem(loads) unfolded_tensor = mod.unfolded_tensor # select lagged tensor and original tensor bundle orig = _(tensor = A.tensor[P+1:T], time_id = K ) # consider the time-dimension of the estimated tensor scalar T = T-P # store estimates bundle Ahat = mod.estimates if exists(meanss) Ahat = tensor_sum(Ahat, meanss) endif # store residuals bundle E = tensor_subtraction(orig, Ahat) if exists(mod.covs) covs = mod.covs else # if LSE was adopted, we use loglik_covs to compute gaussian loglikelihood loglik_covs = whiten_covs(E) endif scalar llik = exists(covs)? Loglik(E, covs) : Loglik(E, loglik_covs) # compute information criteria scalar k=dims[-K]*dims[-K]' scalar tmp=-2*llik scalar aic=tmp+2*k scalar bic=tmp+k*log(T) scalar hqc=tmp+2*k*log(log(T)) # store covariance of residuals matrix e = unfold(E, K) matrix sigma = (1/(T))*e'e # store variance of estimated paremeters: if covs doesn't exist, least squares method is assumed if !exists(covs) csi = compute_variance(&A, loads, sigma, unfolded_tensor) else csi = compute_variance(&A, loads,, unfolded_tensor, covs) endif # Check for stationarity: scalar radius = check_stationarity(loads, dims[-K]) # store standard errors sds = compute_sd_params(csi, dims[-K], P ) # add dimensions matrix est_dims = dims est_dims[K] = T orig.dims = est_dims Ahat.dims = est_dims # add back names of dimensions if they were present if exists(A.dimnames) dimnames = A.dimnames dimnames[K] = dimnames[K][-seq(1,P)] orig.dimnames = dimnames Ahat.dimnames = dimnames endif # store the bundle if !exists(covs) return _(data = orig, estimates = Ahat, param = loads, residuals = E, Sigma=sigma, VCV=csi, sd_param = sds, BIC=bic, loglik = llik, method="LSE", n_lags=P, AIC=aic, HQC = hqc, rho = radius) else return _(data = orig, estimates = Ahat, param = loads, covs = covs, residuals = E, Sigma=sigma, VCV=csi, sd_param = sds, BIC=bic, loglik = llik, method="MLE", n_lags=P,AIC=aic, HQC = hqc, rho = radius) endif /* -doc~ Core function to compute permutation matrix. */ matrix P = zeros(n*m, n*m) loop i = 1 .. n loop j = 1 .. m U = zeros(n,m) U[i,j] = 1 P += U**U' endloop endloop return P /* -doc~ Function to compute the Permutation matrix for k-th dimension */ scalar N = cols(dims) # special case for the first dimension if k == 1 matrix ret = I(prodr(dims)) else matrix P = compute_P(dims[k], prodr(dims[seq(k-1,1)])) endif # special case for the last dimension if k == N ret = P # General case for 1 < k < N elif k != 1 scalar N = cols(dims) scalar dim_iden = prodr(dims[seq(k+1,N)]) matrix ident = I(dim_iden) ret = ident**P endif return ret /* -doc~ Function to compute permutation matrices in all dimensions */ scalar N = cols(dims) matrices Qs = array(N) loop k = 1 .. N Qs[k] = compute_Q(dims,k) endloop return Qs /* -doc~ Function to compute kronecker product of all loadings (or covs) expecpt for the i-th one */ scalar P = nelem(loads) scalar N = nelem(loads[1]) arrays ret = array(P) loop p = 1 .. P matrices ret[p] = array(N) if N==2 ret[p][1] = loads[p][2] ret[p][2] = loads[p][1] else loop i = 1 .. N loads_no_d = loads[p][-i] matrix Phi = {1} loop j = nelem(loads_no_d) .. 1 --decr Phi = Phi**loads_no_d[j] endloop ret[p][i] = Phi endloop endif endloop return ret /* -doc~ This function computes gamma: will be summed to the expected value of the stacked jacobian (Ws) to compute the variance of the estimator */ scalar P = nelem(loads) scalar N = nelem(loads[1]) matrices ret = array(P*(N-1)) # sum of elements of dims squred scalar incr = dims*dims' matrix start_zeros = {} matrix end_zeros = zeros((P-1)*incr,1) scalar count = 0 loop p = 1 .. P if p>1 end_zeros = zeros((P-p)*incr,1) endif # compute gamma for all dimensions except for the last one loop i = 1 .. N-1 count++ matrix a = vec(loads[p][i]) if i ==1 tmp = a | zeros( sum( (dims[-i]).^2 ) , 1) else matrix to_sum_one = dims[1:i-1] matrix to_sum_two = dims[i+1:N] tmp = zeros(sum( (to_sum_one).^2 ),1) | a | zeros(sum( (to_sum_two) .^2) , 1 ) endif ret[ count ] = start_zeros | tmp | end_zeros endloop start_zeros = zeros(p*incr,1) endloop return ret /* -doc~ Function to compute stacked jacobian (Ws in the paper) for all time periods */ scalar N = cols(dims) scalar P = nelem(loads) scalar T = dims[N] matrices Ws = array(T-P) arrays lambda_prod = compute_Phi(loads) loop t = P+1 .. T matrix W_t = {} loop p = 1 .. P matrices Phis = lambda_prod[p][1:N-1] loop i = 1 .. N-1 matrix tmp = unfolded_tensor[t-p][i]*Phis[i]' W_t |= ( tmp **I(dims[i]) ) *Qs[i] endloop endloop Ws[t-P] = W_t endloop return Ws /* -doc~ function to compute H: that is, expected value of the stacked jacobian plus gamma */ scalar T = nelem(Ws) scalar R = rows(Ws[1]) scalar N = nelem(gammas) matrix ret = zeros(R,R) # sigma doesn't exists, LSE covariance if !exists(sigma_inv) loop t = 1 .. T ret += Ws[t]*Ws[t]' endloop else # if sigma does exist, MLE covariance loop t = 1 .. T #ret += Ws[t]*sigma_inv*Ws[t]' ret += qform(Ws[t], sigma_inv) endloop endif # add the adjusted gamma matrix sum_gamma = zeros(R,R) loop i = 1 .. N sum_gamma += gammas[i]*gammas[i]' endloop ret = ret/T return ret + sum_gamma /* -doc~ This functions does the last step in order to compute the variance of the estimator: computes expected value of the stacked jacobian and adds gamma. */ scalar T = nelem(Ws) scalar R = rows(Ws[1]) matrix sum_csi = zeros(R,R) loop t = 1 .. T #sum_csi += Ws[t]*sigma*Ws[t]' sum_csi += qform( Ws[t], sigma) endloop sum_csi = 1/(T)*sum_csi scalar eps = 10e^-9 matrix Hinv = ginv(H, eps ) return (1/T)*(Hinv*sum_csi*Hinv) /* -doc~ Function that handles the pipeline for computation of the variance of the estimator. (see pag. 41 from xi and liao,2021 "Multi-linear Tensor Autoregressive Models") */ matrix dims = exists(A.dims) ? A.dims : get_tensor_dims(A) scalar N = cols(dims) scalar P = nelem(loads) scalar time_id = A.time_id matrices Qs = compute_Qs(dims[-time_id]) matrices Ws = all_Ws(unfolded_tensor, loads, Qs, dims) gammas = compute_gamma(dims[-time_id], loads) if exists(covs) matrix constr_sig = {1} loop i = nelem(covs) .. 1 --decr matrix constr_sig = constr_sig**ginv(covs[i], 10e^-9) endloop matrix H = compute_H(Ws, gammas, constr_sig) ret = compute_csi( H, Ws, constr_sig) else matrix H = compute_H(Ws, gammas) ret = compute_csi( H, Ws, sigma) endif return ret /* -doc~ The function reorders diagonal values of the covariance matrix of the estimator in a set of N matrices of dimensions dims[i] x dims[i], i = 1 .. N, Where N is the number of dimensions excluding time */ scalar K = cols(dims) matrices ret = array(K) scalar ini = 1 scalar fin = 0 loop i = 1 .. K fin += dims[i]^2 ret[i] = mshape(diag_csi[ini:fin],dims[i], dims[i] ) ini = fin +1 endloop return ret /* -doc~ The function applies compute_sd() function to standard errors of estimates of each lag. */ scalar r = rows(csi) matrices variances = msplitby(sqrt(diag(csi)),r/P) arrays ret = array(P) loop p = 1 .. P ret[p] = compute_sd(variances[p], dims) endloop return ret /* -doc~ The function builds label for parameter-matrix of lag p, dimension mat_idx. */ string basic_label = sprintf("A[%d][%d]",p, mat_idx) scalar r = rows(A) scalar c = cols(A) strings labells = array(0) loop j = 1 .. c loop i = 1 .. r labells += sprintf("%s[%d,%d]", basic_label, i, j) endloop endloop return labells /* -doc~ The function prints coefficient, standard error and t-statistic for each value of parameter matrices. We follow vec() order. */ loads = ret.param std_err = ret.sd_param scalar P = nelem(loads) scalar N = nelem(loads[1]) string output = "" loop p = 1 .. P string lag = sprintf("\n\n LAG N. %d:\n", p) output ~= lag loop i = 1 .. N matrix res = vec(loads[p][i]) ~ vec(std_err[p][i]) tmp_mat = loads[p][i] strings names = build_mat_label(tmp_mat, i, p) string dim = sprintf("\n\nMatrix A[%d][%d]:\n",p, i) string output_i = "" outfile --buffer=output_i modprint res names end outfile output ~= dim ~ output_i endloop endloop print output printf "%-18s %-12.7g %-18s %-12.7g\n", "Log-likelihood", ret.loglik, "Akaike criterion", ret.AIC printf "%-18s %-12.7g %-18s %-12.7g\n\n", "Schwarz criterion", ret.BIC, "Hannan-Quinn", ret.HQC print "Check for time series stationarity:" printf "Max. eigenvalue of the Companion Matrix is: %.5f\n",ret.rho if ret.rho>=1 printf "WARNING! The index is more than one, the estimated tenar(%d) model is not stationary and causal.\n",P else printf "The estimated tenar(%d) model is stationary and causal.\n",P endif set verbose off clear --all include TenTS.gfn # Monthly import/export data (in milions of €) for 4 EU nations across 5 sectors, # broken down by EU vs. Extra-EU trade partners. open trade_data.gdt --frompkg=TenTS matrix X = {dataset} # Define dimensions of the tensor scalar npartners = 2 scalar nmonths = $nobs scalar ncountries = 4 scalar nsectors = 5 scalar flow_direction = 2 # Define strings arrays to store dimension names strings time_names = obslabel({obs}) strings countries = defarray("Germany","Spain","France","Italy") strings sectors = defarray("Food", "Crude_mat","Fuels", "Chemicals","Other_man") strings flows = defarray("Imports", "Exports") strings partners = defarray("Extra_EU","EU") # Row-vector containing tensor dimensions matrix dims = {nmonths, ncountries, nsectors, flow_direction, npartners} # Store dimension names in a single array of arrays arrays names = array(cols(dims)) names[1] = time_names names[2] = countries names[3] = sectors names[4] = flows names[5] = partners # Create the tensor: time dimension is the first one scalar time_id = 1 bundle A = tensorize(X, dims, time_id, names) # Sum values across partners (Extra_EU + EU): tensor has now 4 dimensions bundle A = tensor_aggr(A,5,1) # Set time as last dimension (Condition required for tenar_est() ) A = permute(A,{2,3,4,1}) # Compute annual percentage change bundle sdiff_A = tensor_elementratio(tensor_diff(A,12),tensor_lag(A,12)) # Select n. of Lags scalar P = 1 # Estimate TenAr(P) via LSE bundle mod = tenar_est( sdiff_A, P ) # Transform estimated values in levels # if P=1 Alag is from Feb 2002 to Dec 2024 bundle Alag = tensor_drill(A, defbundle("4",seq(P+1,nmonths-12) )) bundle mod.estimates = tensor_sum( tensor_elementprod(mod.estimates,Alag), Alag) # Select observed values that correspond to the fitted ones: # if P=1, A is from Feb 2003 to Dec 2025 A = tensor_drill(A, defbundle("4", seq(12+P+1,nmonths))) # Store data and estimates in a single tensor matrix new_dims = A.dims~2 dimnames = A.dimnames + defarray(defarray("","hat")) bundle res = tensorize(vectorize_tensor(A) | vectorize_tensor(mod.estimates), new_dims, A.time_id, dimnames) # Extract Export and Import bundle export = tensor_drill(res, defbundle("3",2)) bundle import = tensor_drill(res, defbundle("3",1)) # Compute trade balance for both real and fitted data balance = tensor_subtraction(export,import) # PLOT TRADE BALANCE vs ESTIMATES FOR THE 4 COUNTRIES: # Plot observed and fitted trade balance for each country summing values for sectors bundle filt = empty gpbuild balance_plot loop i = 1 .. ncountries filt["1"] = i tensor_plot(balance, balance.time_id ,filt,{2,3},, 0 ) endloop end gpbuild gridplot balance_plot --rows=4 --fontsize=7 --output=display --title="Trade Balance: Real vs Fitted Values"